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ABSTRACT 

Entering a new era of high-energy y-ray experiments, there is an exciting quest for the first 
detection of y-ray emission from clusters of galaxies. To complement these observational ef- 
forts, we use high-resolution simulations of a broad sample of galaxy clusters, and follow 
self-consistent cosmic ray (CR) physics using an improved spectral description. We study CR 
proton spectra as well as the different contributions of the pion decay and inverse Compton 
emission to the total flux and present spectral index maps. We find a universal spectrum of 
the CR component in clusters with surprisingly little scatter across our cluster sample. When 
CR diffusion is neglected, the spatial CR distribution also shows approximate universality; 
it depends however on the cluster mass. This enables us to derive a semi-analytic model for 
both, the distribution of CRs as well as the pion-decay y-ray emission and the secondary 
radio emission that results from hadronic CR interactions with ambient gas protons. In addi- 
tion, we provide an analytic framework for the inverse Compton emission that is produced 
by shock-accelerated CR electrons and valid in the full y-ray energy range. Combining the 
complete sample of the brightest X-ray clusters observed by ROSAT with our y-ray scaling 
relations, we identify the brightest clusters for the y-ray space telescope Fermi and current 
imaging air Cere nkov telescopes (MAGIC, HESS, VERITAS). We reproduce the result in 
Pfrommer d2008l) . but provide somewhat more conservative predictions for the fluxes in the 
energy regimes of Fermi and imaging air Cerenkov telescopes when accounting for the bias 
of 'artificial galaxies' in cosmological simulations. We find that it will be challenging to de- 
tect cluster y-ray emission with Fermi after the second year but this mission has the potential 
of constraining interesting values of the shock acceleration efficiency after several years of 
surveying. Comparing the predicted emission from our semi-analytic model to that obtained 
by means of our scaling relations, we find that the y-ray scaling relations underpredict, by up 
to an order of magnitude, the flux from cool core clusters. 

Key words: magnetic fields, cosmic rays, radiation mechanisms: non-thermal, elementary 
particles, galaxies: cluster: general, Galaxy: fundamental parameters 



1 INTRODUCTION 
1.1 General background 

In the cold dark matter (CDM) universe, large scale structure grows 
hierarchically through merging and accretion of smaller systems 
into larger ones, and clusters are the latest and most massive objects 
that had time to virialise. This process leads to collisionless shocks 
propagating through the intra-cluster medium (ICM), accel erating 
both protons and electrons to highly relativisti c energies dDrurvl 
1 19831 : iBlandford & Eichled[l987l: ISarazin|[l999T) . High resolution 
X-ray observations by the Chandra and XMM-Newton satellites 
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confirmed this picture, with most clusters displaying evidence 
for si gnificant sub s tructur e s, shocks, a n d contact discontinuities 
(e.g. iRosati et all |2002| ; FVoitl 120051 : iMarkevitch & Vikhlininl 
|2007|) . In addition, observations of radio halos and radio relics 
demonstrate the presence of synchrotron emitting e lectrons with 
energies reaching ~ 10 GeV in more than 50 clusters dFerettil2003l : 
iFerrari et al. 2008), although their precise origin in radio halos is 
still unclear. Similar populations of electrons may radiate y-rays 
efficiently via inverse Compton (IC) upscattering of the cosmic mi- 
crowave background photons giving rise t o a fraction of the diffuse 
y-ray background observed by E GRET dLoeb & WaxmarJl200d. 
Totani & Kitavam j200d : lMiniatil2002Ll2003l : llnoue & Nagashimal 
2005h . Although there is no clear observational evidence yet for a 
relativistic proton population in clusters of galaxies, these objects 
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are expected to contain significant populations of relativistic pro- 
tons originating from different sources, such as structure formation 
shocks, radio galaxies, and supernovae driven galactic winds. 



collisions of relativistic protons leading to y-rays ( Volk et al. 


1996; 


EnGlin et al. 1997; Miniati 2003; 


Pfrommer & EnBlin 2003, 


2004; 


Pfrommer et al.l 20081: |Pfrommedl2008l; iKushnir & Waxmad|2009h 


as well as secondary electron injection (Dennison 1980; Vestrand 



Miniati et al] l2001bl: IPfrommer & Enfflml |2004 IPfrommer et all 
20081 ; IKushnir & WaxmarJ 120091 ; IKushnir et all l2009h - The"se 



hadronic collision processes should illuminate the presence of 
these elusive particles through pion production and successive 
decay into the following channels: 

tt* fi* + vJVf, -> e* + vjv e + v M + i>,, 

7T° 2y 

This reaction can only unveil those cosmic ray protons (CRs) 
which have a total energy that exceeds the kinematic threshold 
of the reaction of E,^ = 1.22 GeV. The magnetic fields play 
another crucial role by confining non-thermal protons within the 
cluster volume for longer than a Hubble time, i.e. any protons 
injec ted into the ICM accumulates throughout the cluster's his - 
tory JVolk et alj[l996l ; lEnfilin et alJI 19971 ; llerezinskv et al.lll997h - 
Hence, CRs can diffuse away from the production site, establish- 
ing a smooth distribution throughout the entire ICM which serves 
as efficient energy reservoir for these non-gravitational processes 
jBlasi & Colafrancesco|[l999l : iDolag & Enflliiill2000l : iMiniati et al] 
2001a). 

There is only little known theoretically about the spectral 
shape of the CR population in the ICM. It is an interesting ques- 
tion whether it correlates with injection processes or is signif- 
icantly modified by transport and re-acceleration processes of 
CRs through interactions with magneto-hydrodynamic (MHD) 
waves. The most important processes shaping the CR spectrum 
as a function of cluster radius are (1) acce l eration by structure 
formation shock wav es dOuilis et al. 1 1 19981 ; Miniati et al. 2000; 
IPfrommer et al. 2006t). MHD turbulence, supernova driven galac- 
tic winds dAcciari et al. 2009). or active galactic nuclei (AGN), 
(2) adiabatic and non-adiabatic transport processes, in particular 
anisotropic diffusion, and (3) loss processes such as CR fhermal- 
ization by Coulomb interactions with ambient electrons and catas- 
trophic losses by hadronic interactions. The spectral distribution of 
CRs that are accelerated at structure formation shocks should be 
largely described by a power-law with a spectral index of the one- 
dimensional distribution given by 



r c + 2 
r c -l' 



(1) 



where r c is the shock compression factor. Strong (high Mach 
number) shocks that inject a hard CR population occur either 
at high redshift during the formation of the proto-clusters or to- 
day at the boundary where matter collapses from voids onto fil- 
aments or super-cluster regions. In contrast, merger shocks show 
weak to intermediate str ength with typica l Mach numbers in the 
range of M ^ 2 .4 dRvu et al.l [20031 : IPfrommer et al.l 120061 : 
Skill man et ail 2008). AGNs or supernova remnants a re expected 
to inject CRs with rather flat spectra, a mi as 2.2 - 2.4 jVolk et al.l 
1 1 9961: ISchlickeisej|2002l ; lEnfflirJl2003l) . but it is not clear whether 
they are able to build up a homogeneous population of significant 
strength. 

The CRs offer a unique window to probe the process of struc- 



ture formation due to its long cooling times. While the thermal 
plasma quickly dissipates and erases the information about its past 
history, the CR distribution keeps the fossil record of violent struc- 
ture formation which manifests itself through the spectrum that 
is shaped by acceleration and transport processes. The cluster y- 
ray emission is crucial in this respect as it potentially provides 
the unique and unambiguous evidence of a CR population in clus- 
ters through observing the pion bump in the y-ray spectrum. This 
knowledge enables determining the CR pressure and whether sec- 
ondary electrons could contribute to the radio halo emission. In the 
y-ray regime, there are two main observables, the morphological 
appearance of the emission and the spectrum as a function of po- 
sition relative to the cluster center. The morphology of the pion 
induced y-ray emission should fol low that seen in therm al X-rays 
albeit with a slightly larger extent ( Pfro mmer etai]|2008l) . The pri- 
mary electrons that are accelerated directly at the structure for- 
mation shocks should be visible as an irregular sha ped IC mor- 
phology, most pronou nced in the cluster periphery (Miniatil l2003l : 
IPfrommer e"tai]|2008l) . 



1.2 The y-ray spectrum of a galaxy cluster 

How do the spectral electron and proton distributions map onto the 
y-ray spectrum? We show the CR spectrum within the virial ra- 
dius of a simulated Coma-like galaxy cluster in the upper part of 
Fig.Q] It is shaped by diffusive shock acceleration at structure for- 
mation shocks, adiabatic transport and the relevant CR loss pro- 
cessed- Three distinct features are visible in the spectrum: a cutoff 
close to the proton rest mass at m p c 2 ^ 1 GeV, a concave shape 
for proton energies above m p c 2 and a steepening due to diffusive 
losses at energies E p > 10 16 eV x [/eo/(10 29 cm 2 s~')]~ 3 , where kq 
is the value of the diffusion coefficient at 1 GeV. The dotted lines 
represents different values of the diffusion coefficient which is var- 
ied by a factor two from its fiducial value. The low energy cutoff 
is due to a bala nce of Coulomb a nd hadronic losses at energies 
around a GeV ( Jubel gas et al . 2008). As shown in the present paper 
and in an upcoming work by Pinzke & Pfrommer (in prep.), the 
concave curvature is a unique shap^H that is caused by the cosmic 
Mach number distribution in combination with adiabatic transport 
processes. These features are mapped onto the pion decay y-ray 
emission spectra as a consequence of hadronic CR interactions. 

This can be seen in the lower part of Fig.[T] where the arrows 
indicate the spectral mapping from the CR spectrum to the photon 
spectrum. In a hadronic interaction, CRs produce pions that decay 
into photons with an energy that is on average smaller by a fac- 
tor eight compare to the original CR energy (see Section |2~!2t . At 
CR energies that are larger than the hadronic reaction threshold, 
the CR power-law behavior is linearly mapped onto the pion de- 
cay induced y-ray spectrum (solid blue). This emission component 
clearly dominates the total photon spectrum and therefore shapes 
the total emission characteristics in the central parts of the cluster, 



1 The physics will be thoroughly developed in this work but we will review 
the main characteristics here for introduction. 

2 A CR distribution with uniform spectr al index was also hinted at by 
IMiniati et alj feOOlal) and iMiniatil j2002l) . where the dominating strong 
shocks caused a constant spectral index of about a ^ 2. Their under- 
lying Mach number distribution, however, is in conflict with that ob- 
tained by ind ependent other works j5yu et al. 2003 ; Pfromm er et alj|2007l : 
ISkillman et al.ll2003) and the spectral shape that we find (see Section 18.11 
for more details). 
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Figure 1. Upper Panel: CR spectral distribution within the virial radius of 
a Coma-like cluster. It shows three distinct features: a cutoff around GeV 
energies, a concave shape, and a steepening at high energies due to diffu- 
sive losses. Lower Panel: we show the intrinsic y-ray number flux weighted 
by the photon energy that does not take into account photon propagation 
effects. The arrows indicate the spectral mapping from the CR spectrum 
to the photon spectrum. The pion decay flux is denoted in blue color, the 
secondary inverse Compton (sIC) in red, and the primary inverse Comp- 
ton (pIC) emission in green. Due to the large uncertainty in the diffusion 
coefficient k, we demonstrate how varying k by a factor of two changes 
the corresponding y-ray spectrum at high energies (dotted lines). The green 
band shows how the pIC emission changes if we vary the maximum electron 
injection efficiency from 0.05 (top) to 0.01 (bottom). 



where the densities are high. Note that this spectrum is an intrin- 
sic spectrum emitted at the cluster position and converted to a flux 
while assuming a distance of 100 Mpc without taking into account 
photon propagation effects. Depending on the cluster redshift, the 
finite mean free path of high-energy y-ray s to e + e~-pair produc- 
tion on infra red (IR) and optical photons limits the observable part 
of the spectrum to energies E y <; 10 TeV for clusters with red- 
shifts z ~ 0.03, and sm aller energies for higher redshift objects 
dFranceschini et al]|2008h . 

Secondary CR electrons and positrons up-scatter cosmic mi- 
crowave background (CMB) photons through the IC process into 
the y-ray regime, the so-called secondary inverse Compton emis- 
sion (sIC). This emission component originates from the flat high- 
energy part of the CR spectrum and produces a rather flat sIC spec- 
trum up to the Klein-Nishina (KN) regime. At large electron en- 
ergies, we enter the KN regime of IC scattering where the elec- 
tron recoil effect has to be taken into account. It implies less effi- 
cient energy transfer in such an elastic scattering event compared 
to the Thomson regime and leads to a dramatic steepening of the 
sIC spectrum at y-ray energies around 100 TeV (solid red line). 



The dash-dotted red line shows the hypothetical sIC spectrum in 
the absence of the KN effect (which is never realized in Nature). 
However, it clearly shows that the diffusive CR break is not ob- 
servable in the sIC component for large clusters (while it can move 
to energies below the KN break for small enough clusters, caus- 
ing a faster steepening there). The spectrum shown in green color 
represents the energy weighted photon spectrum resulting from the 
IC process due to electrons accelerated at structure formation and 
merger shocks, the primary inverse Compton emission (pIC). The 
exponential cutoff is due to synchrotron and IC losses which lead 
to a maximum energy of the shock-accelerated electrons. The green 
pIC band shows the effect of the maximum electron injection effi- 
ciency, where we u se an optimistic value of f e ,max = 0.05 (see e.g. 
iKeshet et al]|2003l) in the top and a value of f e ,max = 0.01 at the 
bottom. This more realistic value is suggested to be the theoreti- 
cally allowed upper limit for the injection efficiency that is consis- 
tent with the non-thermal radiat ion of young supernova remnants 
dZirakashvili & AharoniarfeoiOl) . 

This work studies the spectral and morphological emission 
characteristics of the different CR populations in the y-ray regime. 
We concentrate on observationally motivated high-energy y-ray 
bands. (1) The energy regime accessible to the Fermi y-ray space 
telescope with a particular focus on E y = lOOMeV and (2) the en- 
ergy regime accessible to imaging air Cerenkov telescopes (IACTs) 
assuming a lower energy limit of 100 GeV. In Section [2] we de- 
scribe the setup of our simulations, explain our methodology and 
relevant radiative processes considered in this work. In Section [3] 
we study emission profiles and maps, as well as spectral index 
maps. We then present the CR spectrum and spatial distribution 
and show its universality across our simulated cluster sample in 
Section [4] This allows us to derive a semi-analytic framework for 
the cluster y-ray emission in Section [5] which we demonstrate on 
the Perseus and Coma galaxy clusters. Furthermore, we study the 
mass-to-luminosity scaling relations (Section[6j and predict the y- 
ray flux from a large sample of galaxy clusters for the GeV and 
TeV energy regimes in Section [7] We compare our work to previ- 
ous papers in this field and point out limitations of our approach 
in Section [8] We conclude our findings in Section [9] Throughout 
this work we use a Hubble constant of H = 70kms~ 1 Mpc~', 
which is a compromi se between the value found by the Hubble 
key project (Hq = 72. iFreedman et al"] [200 1 ) and from that one in- 
ferred from baryoni c acoustic oscillation measurements (H = 68, 
IPercival et ai1l2009l) . 



2 SETUP AND FORMALISM 

We follow the CR proton pressure dynamically in our simulations 
while taking into account all relevant CR injection and loss terms 
in the ICM, except for a possible proton production from AGN and 
supernova remnants. In contrast, we model the CR electron popu- 
lation in a post-processing step because it does not modify the hy- 
drodynamics owing to its negligible pressure contribution. We use 
a novel CR formalism that allows us to study the spectral properties 
of the CR population more accurately. 

2.1 Adopted cosmology and cluster sample 

The simulations were performed in a ACDM universe using 
the cosmological parameters: Q. m = £1 D m + £\ = 0.3, Ob = 
0.039, fl A = 0.7, h = 0.7, n s = 1, and cr g = 0.9. Here, fl m de- 
notes the total matter density in units of the critical density for 
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geometrical closure today, p clit = 3Hg/(SnG). Ci b , Q.qm and C1 A 
denote the densities of baryons, dark matter, and the cosmological 
constant at the present day. The Hubble constant at the present day 
is parametrized as H = 100 h km s~'Mpc~', while n s denotes the 
spectral index of the primordial power-spectrum, and <xg is the rms 
linear mass fluctuation within a sphere of radius 8 /r'Mpc extrapo- 
lated to z = 0. 

Our simulations were carried out with an updated and ex- 
tended versio n of the distributed-memory parallel TreeSPH code 
GADGET-2 dSpringell 120051 : ISpringel et all 1200 lh . Gravitational 
forces were computed using a combination of particle-mesh and 
tree algorithms. Hydrodynamic forces are computed with a vari- 
ant of the smoothed particle hydrodynamics (SPH) algorithm that 
conserves energ y and entropy where appro priate, i.e. outside of 
shocked regions dSpringel & Hernquistl 2002). Our simulations fol- 
low the radiative cooling of the gas, star formation, supernova feed- 
back, and a photo-io nizing background (details can be found in 
IPfrommer et~aT1l2007l) . 

The clusters have originally b een selected from a low- 
resolution dark-matter-only simulation ( Yoshida et al. 2001). Using 
the 'zoomed initial conditions' technique dKatz & White! 19931) , the 
clusters have been re-simulated with higher mass and force resolu- 
tion by adding short-wavelength modes within the Lagrangian re- 
gions in the initial conditions that will evolve later-on into the struc- 
tures of interest. We analyzed the clusters with a halo-finder based 
on spherical overdensity followed by a merger tree analysis in order 
to get the mass accretion history of the main progenitor. The spheri- 
cal overdensity definition of the virial mass of the cluster is given by 
the material lying within a sphere centered on a local density max- 
imum, whose radial extend R^ is defined by the enclosed threshold 
density condition M(< R a )/(4kR 3 & /3) = p t hres- We chose the thresh- 
old density p C h res (z) = Ap cr i t (z) to be a constant multiple A = 200 
of the critical density of the universe p tT i C (z) = 3H(z) 2 /(8nG). In 
the remaining of the paper, we use the terminology R V1T instead of 
/?2oo- Our sample of simulated galaxy clusters consists of 14 clus- 
ters that span a mass range from 8x 10 13 M o to 3x 10 15 M o where the 
dynamical stages range from relaxed cool core clusters to violent 
merging clusters (cf . Table [T). Each individual cluster is resolved 
by 8 x 10 4 to 4 x 10 6 particles, depending on its final mass. The 
SPH densities were computed from the closest 48 neighbors, with 
a minimum smoothing length (/i sm i) set to half the softening length. 
The Plummer equivalent softening length is 7 kpc in physical units 
after z = 5, implying a minimum gas resolut ion of approximately 

1.1 x 10'°M o (see also lPfrommer et al.ll2007b . 

2.2 Modeling of CR protons and induced radiative processes 

Our simulations follow cosmic ray physics in a self-consistent way 
dPfrommer et alj2006l:lEnfflin et alj2007l ; |jubelgas et alj2008h . We 
model the adiabatic CR transport process such as compression 
and rarefaction, and a number of physical source and sink terms 
which modify the cosmic ray pressure of each CR population sepa- 
rately. The most important source considered^] for acceleration is 
diffusive shock acceleration at cosmological structure formation 
shocks, while the primary sinks are thermalization by Coulomb 
interactions, and catastrophic losses by hadronization. Collision- 
less structure formation shocks are able to accelerate ions and elec- 

3 For simplicity, in this paper we do not take into accou nt CRs injected into 
the inter-stellar medium from supernova remnants (see Aleksic et al. d2010l) 
for a discussion of this topic). 



Table 1. Cluster sample. 



Cluster 


sim.'s 


dyn. state^ 1 ' 


M (2) 

vir 


R l2> 
vir 


kT Q) 

vir 








[M ] 


[Mpc] 


[keV] 


i 
j 


gsa 




Z.O X 1U 


z.v 


1 1 1 


z 


gla 




1 On/ 1 Al5 


Z.J 


1 A £ 


j 


g/za 


rub LLVl 


1 f\ v 10.15 


1 A 
Z.H 


Q A 


A 

4 


g31 


LL 


l.j X 1U 


1 A 
Z.4 


Q A 

y.4 


J 


gib 


Aft 
M 


J.2X 1U 


1. / 


4. / 


6 


g72b 


M 


2.2 x 10 14 


1.2 


2.4 


7 


glc 


M 


2.0 x 10 14 


1.2 


2.3 


8 


g8b 


M 


1.5 x 10' 4 


1.1 


1.9 


9 


gld 


M 


1.3 x 10 14 


1.0 


1.7 


10 


g676 


CC 


1.3 x 10 14 


1.0 


1.7 


11 


g914 


cc 


1.2 x 10 14 


1.0 


1.6 


12 


gle 


M 


9.1 x 10 13 


0.93 


1.3 


13 


g8c 


M 


8.5 x 10 13 


0.91 


1.3 


14 


g8d 


PreM 


7.8 x 10 13 


0.88 


1.2 


Notes: 



(1) The dynamical state has been classified through a combined crite- 
rion invoking a merger tree study and the visual inspection of the X- 
ray brightness maps. The labels for the clusters are M-merger, PostM- 
post merger (slightly elongated X-ray contours, weak cool core (CC) re- 
gion developing), PreM-pre-merger (sub-cluster already within the virial 
radius), CC-cool core cluster with extended cooling region (smooth X- 
ray profile). (2) The virial mass and radius are related by Ma(z) = 
|ffAp cr j[(z)S 3 , where A = 200 denotes a multiple of the critical over- 
density p cr it(z) = 3H(z) 2 K%nG). (3) The virial temperature is defined by 
&7vir = GAfyir// m p /(2i? v ir), where /i denotes the mean molecular weight. 



trons in the high-energy tail of their Maxwellian distributio n func- 
tions through diffusive shock acceleration (for reviews see | Drurvl 
1 19831: iBlandford & Eichlerll 19871 : iMalkov & O'C Drurvll200lh . In 
the test particle picture, this process injects a CR distribution with 
a power-law in momentum and a slope that depends on the instan- 
taneous sonic Mach number of the shock. The overall normaliza- 
tion of the injected CR distribution depends on the adopted sub- 
resolu tion model of diffusive shock acceleration (e.g., EnBl in et"al] 
120071) : in particular it depends on the maximum acceleration ef- 
ficiency fmax,p = ScR.max/fidiss which is the maximum ratio of CR 
energy density that can be injected relative to the total dissipated 
energy density at the shock. We assume that in the saturated regime 
of shock acceleration, 50 percent of the dissipated energy at strong 
shocks is injected into cosmic ray protons. While there are indi- 
cations from super nova remnant observations of one rim region 
llHelderetalJl2009l) as well as theoretical studies dKang & Jonesl 
2005) that support such high efficiencies, to date it is not clear 
whether these efficiencies apply in an average sense to strong colli- 
sionless shocks or whether they are realized for structure formation 
shocks at higher redshifts. This high efficiency rapidly decreases for 
weaker shocks (decreasing Mach n umber) and eventua lly smoothly 
approaches zero for sonic waves dEnfilinetal1l2007h . Our paper 
aims at providing a quantitative prediction of the y-ray flux and 
hence the associated CR flux that we expect in a cluster depend- 
ing on our adopted acceleration model. Non-detection of our pre- 
dicted emission will limit the CR acceleration efficiency and help 
in answering these profound plasma astrophysics questions about 
particle acceleration efficiencies. 

We significantly revised the CR methodology and allow for 
multiple non-thermal cosmic ray populations of every fluid element 
(Pinzke & Pfrommer, in prep.). Each CR population f,{p,R) is a 
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power-law in particle momenturrQ, 
f i (p,R) = C i (R)p- a '6(p-q i ), 



(2) 



characterized by a fixed slope a h a low-momentum cutoff q h and 
an amplitude C,(/J) that is a function of the position of each 
SPH particle through the variable R. For this paper we have cho- 
sen five CR populations with the spectral index distribution a = 
(2.1,2.3,2.5,2.7,2.9) for each fluid element (a convergence study 
on the number of CR populations is presented in the appendix I A2t . 
This approach allows a more accurate spectral descriptioijf] as the 
superposition of power-law spectra enables a concave curvature of 
the composite spectrum in logarithmic representation. Physically, 
more complicated spectral features such as bumps can arise from 
the finite lifetime and length scale of the process of diffusive shock 
acceleration or incomplete confinement of CRs to the acceleration 
region. These effects imprint an upper cutoff to the CR population 
locally that might vary spatially and which translates into a convex 
curvature in projection. Additionally, interactions of pre-existing 
CRs with MHD waves can yield to more complex spectral features. 
Future work will be dedicated to study these topics. 

In addition to the spectral features mentioned above, we model 
in the post-processing the effect of high-energy CR protons that are 
no longer confined to a galaxy cluster as these are able to diffuse 
into the ambient warm-hot intergalactic medium (WHIM). In this 
paper we define WHIM to be the regi on within /? V j r < R < 3 R V!Y , 
which is a subset of the entire WHIM (Dav e et al .1200 ll) . Assuming 
particle scattering off magnetic irregularities with the Kolmogorov 
spectrum, we obtain the characteristic scaling of the diffusion coef- 
ficient k ^ kq (E/Eq)' 13 , where we normalize k at £o = 1 GeV. One 
can estimate the c haracteristic proton energy E p break at wh ich the 
spectrum steepens dVolk et al.fl996l ; lBerezinsky et aljl997b . 



J p, break 



(6kot) j 



10" GeV 



(3) 



\3Mpc) llO^cnr's- 1 / \T Hum J 



For the reminder, we adopt a value of the diffusivity that is scaled to 
R = 2 R vir for each cluster, as this volume is expected to fall within 
the virialised part of the cluster past the accretion shock region 
jPfrommer et all 2007; Mol nar et all2009T ) and traps CRs in a clus- 
ter for time scales longer than a Hubble time. This choice also has 
the property that the diffusion break is at energies E p > 100 GeV; 
hence it does not interfere with the pion decay as well as secondary 
IC emission in the energy regime accessible to IACTs as we will 
show in the following. The momentum of a photon that results from 
pion decay is given by 

K p P p P p 

P y ^— — «— . (4) 

r 2 £ 8 

This approximate relation is derived using the inelasticity K p « 1/2 
and multiplicity f ss 2 for the p + p — > n° channel together with 
the two photons in the final state. Secondary electrons that are in- 
jected in hadronic CR interactions Compton up-scatter CMB pho- 
tons. A break in the parent CR spectrum would imprint itself in 



4 The true CR particle momentum is denoted by P p , but we loosely refer 
to p = P p l(m p c) as the particle momentum. 

5 The total CR proton spectrum is a sum of the spectra of the individual 
SPH particles within a certain volume, and since our sample contains a 
large number of SPH particles with varying normalization, the CRp spectral 
index is a statistically well defined continuous quantity. 



the sIC spectrum if there are no other effects that modify the spec- 
trum at lower energies. Compared to the pion decay emission, this 
break manifests at slightly higher energies (for parameters adopted 
in Fig.|T|l. The momentum of the electrons P c depends on the proton 
momentum P p through the relation given by hadronic physics 



p „ £p£p „ £p 

c ~ 4 ^ ~ 16' 



(5) 



Here we used the p + p — > i& channel together with the four parti- 
cles in the final state of the charged pion decay (e ± , v e /v e , v p , v p ). 
Combining the classical inverse Compton formulae from CR elec- 
trons with energies E e > 1 GeV 



4 F 

-£CMB 



(6) 



with the energy relation in equation (|5j we obtain a break in the 
sIC spectrum. This steepening of the CR spectrum take place at 
high photon energies ,E s ic, break - 10 17 eV where we choose CMB 
photons with the energy £Vjmb = ^ y CMB - 0.66 meV as source 
for the inverse Compton emission using Wien's displacement law. 
It turns out that these energies are deeply in the Klein-Nishina 
regime. This means that in the rest frame of the energetic electron, 
the Lorentz boosted photon energy is comparable to or larger than 
the electron rest mass, Eic = yjiv lv i t ~ m e c 2 , so that the scattering 
event becomes elastic. This implies a less efficient energy trans- 
fer to the photon and manifests itself in a break in the resulting 
IC spectrum. While the number flux scales as f ~ E ] 



-(Be-l)/2 
IC 



the Thomson-limit for £ lc <k 30 TeV, it steepens significantly to 
T ~ £ K f c log(£ IC ) in the extreme KN-limit for E [c » 30 TeV, 
whe re q e is the spectral index of the (cooled) CR electron distribu- 
tion dBlumenthal & Gouldll970|) . 



2.3 Magnetic fields 

High energy CR electrons with y c > 200 loose their energy by 
means of IC scattering off CMB photons as well as through inter- 
actions with cluster magnetic fields which results in synchrotron 
emission. The relative importance of these two emission mecha- 
nisms depends on the rms magnetic field strength, 6, relative to 
the equivalent field strength of the CMB, B C mb = 3.24(1 + zfpG, 
where z denotes the redshift. In the peripheral cluster regions, 
where B <k Scmb. the CR electrons loose virtually all their energy 
by means of IC emission. In the central cluster regions, in particu- 
lar in the dense centers of cool cores, the magnetic energy density 
is probably comparable or even larger than the energy density of 
the CMB JVogt & Enfflinll20051) . e ph = B 2 CMB /(Sn). Hence in these 
regions, the radio synchrotron emission carries away a fraction of 
the CR electrons' energy losses; an effect that reduces the level of 
IC emission. We model the stren gth and morphology of the mag- 
netic fields in the post-processing ( IPfrommer etai] 2008) and scale 
the magnetic energy density field e B by the thermal energy density 
6 th through the relation 



2a B 



(7) 



where e Bo = B?/(8?r) and e lho denote the core values. If not men- 
tioned otherwise, we use the magnetic decline a B = 0.5 and the 
central magnetic field So = 10yuG throughout this paper. The cen- 
tral thermal energy density e lho = 3P t h /2, is calculated by fitting 
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the modified /3-model 



P{R) 



thq 



1 + 



R 

*^core 



(8) 



to the radial pressure P(R). The parametrization in equation (7} 
is motivated by both cosmological MHD SPH simulations 



jPolag et alj 19991) and radiative ada ptive mesh refinement MHD 
simulations jDubois & Tevssiej2 008). Rather than applying a den- 
sities scaling as those simulations suggest, we use a scaling with 
thermal gas energy density which is not affected by the over-cooled 
centers in radiative simulations that do not take into account AGN 
feedback. 



2.4 CR electron acceleration and inverse Compton emission 

2.4.1 Modeling diffusive shock acceleration 

Collisionless cluster shocks are able to accelerate ions and electrons 
through diffusive shock acceleration (for reviews see jDrurvll 19831 ; 
Blandf ord & Eichlerlll987l ; iMalkov & O'C Drury||20oTir Neglect- 
ing non-linear shock acceleration and cosmic ray modified shock 
structure, the process of diffusive shock acceleration uniquely de- 
termines the spectrum of the freshly injected relativistic electron 
population in the post-shock region that cools and finally dimin- 
ishes as a result of loss processes. The y-ray inverse Compton emit- 
ting electron population cools on such a short time scale r sync < 



10 yrs (compared to the long dynamical time scale T iyn ~ 2 Gyr) 
that we can describe this by instantaneous cooling^] In this approx- 
imation, there is no steady-state electron population and we would 
have to convert the energy from the electrons to inverse Compton 
and synchrotron radiation. Instead, we introduce a virtual electron 
population that lives in the SPH-broadened shock volume only; this 
is defined to be the volume where energy dissipation takes place. 
Within this volume, which is co-moving with the shock, we can use 
the steady-state solution for the distribution function of relativistic 
electrons and we assume no relativistic electrons in the post-shock 
volume, where no energy dissipation occurs. Thus, the cooled CR 
electron equilibrium spectrum can be derived from balancing the 
shock injection with the IC/synchrotron cooling: above a GeV it is 
given by 

P 



f e (Pt) = C e P e 



C, oc 



(9) 



Here, we introduced the dimensionless electron momentum p e = 
P e /(m e c), where P e is the electron momentum, a e = a m j + 1 is the 
spectral index of the equilibrium electron spectrum, p is the gas 
density, e b is the magnetic energy density, and £ ph denotes the pho- 
ton energy density, taken to be that of CMB photons. The primary 
CRe distribution in equation l[9} is calculated in the post-processing 
with a spectrum reflecting the current Mach number of the shock 
(without the assumption of spectral bins). Superposing the individ- 
ual spectra of a large number of SPH particles, each with a spec- 
trum reflecting the accelerating shock, produces a well defined total 
spectrum with a running index in general. A more detailed discus- 
sion o f this simplified approach can be found in Pfrommer et al. 
d2008l) . 

Once the radiative cooling time due IC and synchrotron 
emission becomes comparable to the diffusive acceleration time 



6 Assuming a magnetic field of a few fiG and an electron dens ity n e = 
10~ 3 c m -3 , for further discussion see e.g. appendix in Pfrommer et al. 
120081) 



scale, the injection spectrum experiences a high-energy cutoff 
(Web b G.lvi] Il984l) . The electrons start to pile-up at this crit- 
ical energy; the super-exponential term describing the maxi- 
mum energy of electrons reached in this process, however, ef- 
fectively cancels this pile-up feature which results in a pro- 
longed power-law up to the elec tron cutoff momentum p c ~ p m:lx 
dZirakashvili & Aharonianll2007t) . We account for this effect by us- 
ing the following parametrization of the shock injected electron 
spectrum, 



fe(*, Pc) 



CePe" mJ [1 + eXp 



(10) 



where x is the distance from the shock surface, j(x,p c ) and S s de- 
scribe the characteristic momentum and shape of the pile-up region. 
The continuous losses cause the cutoff to move to lower energies 
as the electrons are transported advectively with the flow down- 
stream. Integration over the post-shock volume causes the cutoffs 
to add up to a new power-law that is steeper by unity compared to 
the injection power-law (equation |5J. Hence, the shock-integrated 
distribution function - as defined in equation dB 13b and displayed 
in Fig. IB II - shows three regimes: (1) at low energies (but large 
enough in order not to be affected by Coulomb losses) the original 
injected power-law spectrum, (2) followed by the cooled power-law 
that is steeper by unity, a e = a m j + 1, and (3) an ultimate cutoff that 
is determined from the magnetic field strength and the properties 
of the diffusion of the electron in the shock (we refer the reader to 
Section lB2l for a detailed discussion). Note that only the last two 
regimes are important for y-ray IC emission. 

The fact that we observe X-ray synchrotron emission at shocks 
of young supernova remnants dVink et al.ll2006l ; ISlane et alj|l999l) 
necessarily requires the existence of high-energy CR electrons with 
E c ^ 25 TeV. The non-thermal synchrotron emission generated by 
CR electrons with energies E c > GeV is given by 



/; v. 



synch ' 



3eBh 
2nrruc 



y 



(ii) 



where y ~ 5 x 10 7 (E e ~ 25 TeV) and magnetic fields of order 
100 (jG are required to reach X-ray energies of order 10 keV. To 
keep the highly relativistic electrons from being advected down- 
stream requires efficient diffusion so that they can diffuse back up- 
stream crossing the shock front again. We therefore use the most 
effective diffusion, refereed to as Bohm diffusion limit, as the elec- 
tron propagation model at the shock. Balancing Bohm diffusion 
with synchrotron/IC cooling of electrons enables us to derive a 
maximum energy of the accelerated CR electrons at the position of 
the shock surface (derived in appendix, see also IWebbGJyfll 19841; 
lEnfilin et all 19981) 



Be u m e c T[ 0SS 



1 



rc(r c + l) 



(12) 



where diffusion parallel to the magnetic field has been assumed. 
Here u denotes the flow velocity in the inertial frame of the shock 
(u = v s ), and the electron loss time scale due to synchrotron and 
inverse-Compton losses reads 



4cr T y 

3ra c c 



(s B + e p h) , 



(13) 



where cr T is the Thompson cross-section. r c denotes the shock com- 
pression and is given by 



(y, h + l)M 2 
(y th -l)Af + 2 



(14) 
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Table 2. Electron energy cutoff in different regions of a typical cluster. 



Cluster variables' 1 ' 


< 0.03 R vil 


0.3R vil - 


Rvi! 


WHIM 


/i sm l [kpc] 


< 30 


60 


120 


240 


n[10~ 3 X particles cm -3 ] 


4 


0.4 


0.04 


0.004 


r c 


1.3 


2.3 


3.0 


3.7 


M 


1.2 


2 


3 


6 


T [keV] 


15 


9 


5 


0.2 


c SO und[kms~'] 


2400 


3100 


3500 


1400 


B IjuG] 


10 


2.5 


0.6 


0.04 


£max [TeV] 


50 


100 


65 


6.5 



Notes: 

(1) Other constants used: mean molecular weight /i = 0.588 m p . 



where M denotes the sonic Mach number. Inserting typical num- 
bers for different cluster regions show that the electron cutoff en- 
ergy, £ max , only varies within a factor two inside ^ v j, (Table[2}. The 
equivalent cutoff energy in the IC spectrum can easily be derived 
from equation ((6), yielding Ei Ccut ^ 10 TeV x (£ m!lx /50TeV) 2 . 



2.4.2 IC emission 

Following iRvbicki & Lightmarj Jl979t) . we calculate the inverse 
Compton emission from electrons that up-scatter CMB photons. It 
should be noted that we neglect the inverse Compton emission in- 
duced by starlight and dust, which might contribute significantly in 
the inner 10 kpc of the cluster. The inverse Compton source density 
A\c in units of produced photons per unit time interval and volume 
for a simple power-law spectrum of CRes (equation[9j scales as 



(15) 



where a v = (a e - l)/2. When we account for the competition 
between radiative cooling and diffusive acceleration of electrons 
in the shock region, the shape of A K in the high-energy regime 
changes. Using the cooled electron distribution of equation JB 13b . 
we construct an effective integrated source function for primary in- 
verse Compton emission (see Section |B3l for a self-consistent and 
extensive description), 



tpIC 



^0(<Tc.raax, C e ) /ic(tte) /kn(£iC, <*e) 



kT c 



1+0.84 



£ic,cu( 



\<5lc(£lC.O'c) 



exp 



4.07 E lc 

£ic,cut 



(16) 



where Bohm diffusion has been assumed. The normalization con- 
stants Ao(£ c ,max> C e ) an d /ic(Q'c) are derived in equations (B201 and 
dB 19b . respectively. The KN suppression of the IC spectrum is cap- 
tured by /kn(£ic, a e) (equation IB22b . and the shape of the tran- 
sition region from the Thompson to the KN regime is given by 
<5ic(£ic, <*e) (equation IB24b . Following recent work that carefully 
models the non-thermal radiatio n of young supernova remnants 
dZirakashvili & Aharon ian 2010), we typically adopt a maximum 
electron injection efficiency of f c ,max = 0.01. We note that this value 
seems to be at the upper envelope of theoretically allowed values 
that match the supernova data. 



2.5 Multiphase structure of the ICM 

The ICM is a multiphase medium consisting of a hot phase which 
attained its entropy through structure formation shock waves dis- 
sipating gravitational energy associated with hierarchical cluster- 
ing into thermal energy. The dense, cold phase consists of the 
true interstellar medium (ISM) within galaxies and at the cluster 
center as well as the ram -pressure str ipped ISM that has not yet 
dissociated into the ICM dPolag et alj|2009h . All of these phases 
contribute to the y-ray emission from a cluster. Physically, the 
stripped ISM should dissociate after a time scale that depends on 
many unknowns such as details of magnetic drap ing of ICM fields 
(Dursi & Pfrommer 2008; P frommer & Dursill2009h on galaxies or 
the viscosity of the ICM. In SPH simulations, this dissociation pro- 
cess is suppressed or happens only incompletely in our simulations 
leaving compact galactic-sized point sources that potentially bi- 
ases the total y-ray luminosity high. On the other hand, once these 
stripped compact point sources dissociate, the CRs diffuse out in 
the bulk of the ICM, and produce y-rays by interacting with pro- 
tons of the hot dilute phase. This flux, however, is negligible since 



oc J dV reaiff.-CRs "icm « J dV n 



ism— CRs "ism 



"icm— C 



, where — ~ 10 . (17) 



Here n lcm (n 1!>m ) denotes the gas number density in the ICM (ISM), 
M ism-CRS the CR number density in the ISM, and n diff _ CRS the CR 
number density of the CRs that diffused out of their dense ISM 
environment into the ambient ICM that is in pressure equilibrium 
with the ISM. In the second step of equation (T7J we accounted for 
the adiabatic expansion that these CRs would experience as they 
diffused out. In the last step we assumed Li sm _cR S ~ iicm-CRs, i-e. 
that the CR luminosity of all compact galactic sources in a cluster 
is of the same order as the CR luminosity in the diffuse ICM; a 
property that is at least approximately true in our simulations as we 
will show later on. Leaving all gaseous point sources would defini- 
tively be too optimistic, removing all of them would be too con- 
servative since cluster spiral galaxies should contribute to the total 
y-ray emission (which defines our so-called optimistic and conser- 
vative models). Hence, we perform our analysis with both limit- 
ing cases, bracketing the realistic case. The effect from the gaseous 
point sources is largest in low mass clusters, where they constitute a 
few percent of the total ICM mass. For high-mass clusters this frac- 
tion is lower, and constitutes only about one percent. In practice, we 
cut multiphase particles with an electron fraction x c < 1.153 and a 
gas density above the star forming threshold 2.8 x 10~ 25 gcnT 3 . 
If nothing else is stated, we use our conservative model without 
galaxies throughout the paper. 



3 CHARACTERISTICS OF y-RAY EMISSION 

From surface brightness S y maps that are obtained by line-of-sight 
integration of the source functions, we study the y-ray emission to 
characterize the morphology of clusters. Additionally, we use emis- 
sion profiles to compare pIC, sIC, and pion decay induced emission 
for different clusters. 



3.1 Morphology of y-ray emission 

The left side of Fig.|2]shows the morphology of the y-ray emission 
above 100 GeV that results from hadronic CR interactions with am- 
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Figure 2. The y-ray emission above lOOGeV of our Coma-like cluster g72a is shown. We show the pion decay y-ray emission that originates from 
hadronic CR interactions with ambient gas protons on the left. On the right, we plot the inverse Compton emission from both, primary and secondary 
CR electrons. Comparing the different y-ray emission components, we note that the pion decay has a very regular morphology and clearly dominates 
the cluster region. In contrast, the emission from primary electrons shows an irregular morphology that traces the structure formation shock waves and 
dominates in the virial periphery and the warm-hot intergalactic medium. 



tor* 




5 




Figure 3. Projected energy dependent photon index for the galaxy cluster g72a. On the left, we show the photon index between 100 MeV - 1 GeV and 
on the right, the photon index between 100 GeV - 1 TeV. Note that the central part (where pion decay dominates) shows little variations which reflects 
the spectral regularity of the CR distribution. At the periphery and beyond, there are larger fluctuations due to the inhomogeneity of the pIC emission. 



bient gas protons. The right side shows the primary and secondary 
IC emission for the post-merger cluster g72a. The comparison of 
the two panels shows that the central parts are dominated by the 
pion induced y-ray emission. It has a very regular morphology that 
traces the gas distribution. There is a transition to the pIC emis- 
sion as the dominant emission mechanism outside the cluster in the 
WHIM at a level depending on the dynamical state of the cluster. 
The pIC emission is very inhomogeneous which can be easily un- 
derstood since it derives from primary CR electrons that are directly 
accelerated at structure formation shocks. Structure formation is 
not a steady process, it rather occurs intermittently. The morphol- 
ogy of the y-ray emission above 100 MeV for g72a was investigated 



in lPfrommer et alj d2008h ■ It shows a very similar morphology, in- 
dicating a similar power-law CR spectrum. 

The projected energy dependent photon index is a well defined 
continuous quantity as it is defined through 



dlogS y 
dlog£ r 



lot 



r£(*x) = -- 



S 7 (E h R ± ) 
S 7 (E 2 ,R ± ) 



Ml) 



(18) 



where S y (E y ,R ± ) is the surface brightness of y-rays with energies 
> E y at the projected radius R ± . Using 5 r maps, such as the ones 
in Fig. [2] we can extract the photon index V E 2 between the energies 
Ei to Z?2- Close to the pion bump (see Fig.QJ at m n o/2 ^ 68.5 MeV 
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Figure 4. Azimuthally averaged profiles of the y-ray emission components above 100 MeV for two different clusters: large cool core (CC) g8a (left 
panel) and large merger g72a (right panel). Shown are the dominating pion decay induced emission (dash-dotted lines), the IC emission of primary 
CR electrons accelerated directly at structure formation shocks (dotted lines), the secondary IC electrons resulting from CRp-p interactions (dashed 
lines) and the sum of all emission components (thick line). We cut the galaxy emission from our individual emission components in both panels and 
additionally show the total emission including galaxies (grey solid line). 



(energy of a photon originating from a decaying n° at the thresh- 
old of the hadronic p-p reaction) the photon spectrum has a convex 
shape. This is characterized by a flatter photon index compared to 
the asymptotic limit, E y » m n ac 2 . Calculating r^j at two discrete 
energies results in a slightly steeper value for T^ 2 than its continu- 
ous counterpart r at E\ . 

Since we are interested in comparing the morphology of clus- 
ters to spectra, we calculate the projected photon index for g72a 
(see Fig. [3j. We compare T% for an energy regime accessible to 
the Fermi space telescope of 100 MeV - 1 GeV (left panel) to 
100 GeV - 1 TeV accessible for IACTs (right panel). In the Fermi 
regime, we find a median value of r],^ eV =^ 0.9 in the central re- 
gions of the cluster and a value rj,^ cV ^ 1.1 in the periphery. 
The reason is that the total emission in the central regions of the 
cluster is dominated by the pion decay emission at 100 MeV with 
a lower spectral index than pIC due to the pion bump. In the pe- 
riphery and the WHIM, where the pIC contributes substantially 
to t he total emission, inte rmediate shocks with M ~ 4 are typi- 
cal jPfrommer et al.l2006T) . Using the spectral index of the electron 
equilibrium spectrum, a e = a iu j + 1, where a m j m 2.2 for M — 4 
and Sic ~ £^T e ~ — E^J' 1 , results in the observed steepening of 
r !ooMeV in the periphery. 

We now turn to the energy region important for IACTs with 
the photon index r[^J^ v . In the central regions the photon index 
traces the proton spectral index r|^J^ v = a - 1 ~ 1.25 since this 
spatial region is dominated by the asymptotic regime of the pion 
emission. Moving towards the periphery, the photon index steepens 
to r|^Q cV > 1.4, despite the presence of strong external shocks 
as well as accretion shocks that efficiently accelerate electrons on 
these large scales. The reason for this steepening is the exponential 
cutoff in the pIC emission. Increasing the energy or the distance 
from the cluster results in an even steeper photon index. 



3.2 Emission profiles 

The profiles of different non-thermal y-ray emission processes 
without galaxies are shown in Fig. [4] for, a large CC cluster (g8a, 
left) and large post-merging cluster (g72a, right). The secondary 



IC emission traces the dominating pion decay emission because to 
zeroth order, both components depend on n gils n C R, where n gas is the 
gas number density and n C R the CR number density. This would be 
exactly true if the magnetic field was smaller than Bcmb - 3trG. 
In this case, the CRe population would exclusively cool by means 
of IC emission. Since we assume the central magnetic field to be 
larger than Bcmb, a fraction of the CRe energy is radiated through 
synchrotron emission into the radio, causing a larger discrepancy 
of the sIC emission compared to the pion emission in the center. 

In contrast to the centrally peaked secondary emission compo- 
nents, the average primary IC emission shows a rather flat surface 
brightness profile which can be nicely seen in our cool core clus- 
ter g8a. This is because we see the strong accretion shocks that 
efficiently accelerate CRes (in terms of the injected energy den- 
sity) in projection. There are noticeable exceptions in the centers 
of both clusters: accreting small sub-clumps dissipate their grav- 
itational energy through weak shocks in the larger core regions 
of clusters. However, once these weak shock waves encounter the 
(over-cooled) centers of our simulated clusters they transform into 
strong (high Mach number) shock waves. These inject a hard pop- 
ulation of primary CR electrons which causes the centrally peaked 
and bright pIC emission. We also observe an excess pIC emission 
at a radius of ~ 1 Mpc in g72a. This can be traced back to a promi- 
nent merger shock wave with a Mach number up to M — 3.5 that 
accelerates primary CRes (see also Fig.|2j. 

The total emission flattens out in the cluster exterior close to 
Z? v j r because of two reasons: (1) there the pIC emission contributes 
significantly to the total emission because of strong merger shocks 
as well as accretion shocks, and (2) subhalos in the periphery that 
have not yet merged with the larger halo contribute to the pion 
decay induced emission. In this regime, the halo-halo correlation 
term starts to dominate the average density profile o f a cluster with 
its characteristic flattening lHavashi"& White 2008). This naturally 
translates to the pion emission profile that tightly correlates with 
the gas density distribution. 

The ratio between pion decay and sIC emission can be esti- 
mated analytically by calorimetric considerations. To this end we 
compare the CR energy spectrum, E 2 f(E), at the respective CR 
energies which give rise to y-ray emission at some specific pho- 
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Figure 5. The y-ray number flux weighted by photon energy of our Coma-like cluster g72a: core region (left panel) and WHIM region (right panel). The 
pion decay flux is shown in blue color, the secondary IC in red, and the primary IC emission in green. The solid IC lines assume a central magnetic field 
strength of Bo = IOjuG while the dashed lines assume Bq = I fiG, The lower panels show the photon spectral index T for each emission component, 
where the colors are the same as in the upper panel. Note that in the core region (R < 0.3i? v j r ) the pion decay induced flux dominates the emission. In 
the WHIM region (R v i r < R < 3 i? v j, ) the pion decay and sIC emission components are much smaller compared to the central part. Contrary, the pIC 
component is boosted in the WHIM due to accretion shocks onto our simulated cluster. Nevertheless, its emission falls slightly short of the pion decay. 



ton energy, say E y = 1 GeV. Additionally we have to include a 
factor, / s , that accounts for the possibility that the CRes do not 
only emit IC y-rays but also radio synchrotron radiation. Using the 
hadronic branching ratios for the production of pions, Rj> ^1/3 
and R„± - 2/3, as well as the multiplicities for the decay products 
in the respective decay channels, £ n o^ y — 2 and - 1/4, we 

obtain 



E 1 s 



^p.sIC sIC 



R 7fi £r°- 
£r±- 



-p.sIC 



/b - 30 / B . 



(19) 



Here <r^ c = 2.3 is the CR proton spectral index between the CR 
energy E vn o = 8.0 GeV that give rise to pion decay flux at 1 
GeV and the CR energy £ p . s ic = 8.0 TeV that gives rise to sIC 
flux at 1 GeV. The y-ray so urce function for pion decay and sIC 
rtPfrommer & Enfflinl |2004)) is denoted by s^> and s ic , respec- 
tively. Finally, the factor f B = (B/Bcmb) 2 + 1 — 1, for magnetic 
fields much smaller than the CMB equivalent magnetic field, oth- 
erwise /b > 1. In the region close to 03R YII , the magnetic field is 
about 2.4 pG in our model (Table [2}, which implies an emission 
ratio of about 50. For smaller photon energies than 1 GeV, the pion 
decay y-ray emission falls below the asymptotic power-law due to 
the characteristic pion bump. This effect implies a lower ratio of 
about 40 (instead of the expected ratio of about 100) for the emis- 
sion above 100 MeV in Fig. [4] 

We also show the total surface brightness profile with galax- 
ies in Fig. [4] (grey line). The resulting profile shows a boosted 
emission by about a factor two compared to the one where we ex- 
clude galaxies from the surface brightness. The entire population of 
these galaxies takes up only a negligible volume so that the volume 
weighted CR pressure is almost the same in either case, when tak- 
ing these galaxies into account or not. We note that this bias needs 
to be addressed when deriving av erage CR pressure c ontributions 
from the cluster's y-ray emission jAleksig et al]|2O10h , Especially 
in the inner parts, the profile is very inhomogeneous. Since galax- 
ies follow an approximate Poisson distribution and since the inner 
radial bins of the profile sample only few galaxies, we naturally 
obtain a larger Poissonian scatter across the inner radial bins. 



3.3 Emission spectra from the cluster core and WHIM 

The central parts of clusters are characterized by high gas and CR 
densities, and magnetic fields - at least compared to average val- 
ues of the ICM. Even though the cluster core region only makes 
up a fraction of the total volume of the ICM, the high densities re- 
sult in a significant y-ray flux contribution to the total flux from 
the cluster. In contrast to the cluster center, the WHIM is character- 
ized by on average low gas and CR densities, and magnetic fields. 
The low densities cause a smaller total y-ray flux from this region 
compared to the cluster core regions. However, the WHIM of the 
super-cluster region contains a large number of galaxies and groups 
that are accreted onto the cluster. This generates more shocks com- 
pared to the cluster core region. The cluster characteristics in the 
two regions give rise to different normalizations of the individual 
y-ray emission components, but with a similar shape. The shape 
of the emission components from the different regions agrees with 
that of the entire cluster as shown in Fig. [T] The emission can be 
summarized as follows. The ;r°-decay is characterized by the so- 
called pion bump followed by a concave curvature and a diffusive 
break. The pIC emission component shows a power-law followed 
by an exponential cutoff while the sIC component has a power-law 
with similar index that is however followed by the Klein-Nishina 
break. 

In Fig. [5] we show the y-ray number flux weighted by pho- 
ton energy from different regions of our Coma-like cluster g72a 
that we place at the distance of 100 Mpc. The left figure shows 
the y-ray number flux within the core region (if < 0.3 R vu ), where 
the n° '-decay dominates over the sIC component that itself is sub- 
dominant to the pIC component. The surface brightness profile of 
the 7r°-decay is sufficiently flat in the core region so that the y-ray 
flux is dominated by the outer scale around if ~ 0.3 if v i r where most 
of the volume is. Hence the 7r°-decay flux is largely insensitive to 
numerical inaccuracies of our modeling of the physics at the very 
center of the cluster. Both the pIC and sIC emission components 
have a larger y-ray flux in our models with weak central magnetic 
fields (B = 1 pG, /b - 1) compared to our models with strong cen- 
tral magnetic fields (Bq = 10//G, /b — 2). The sIC emission with 
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S = lOytiG is characterized by B =s B C mb on scales around the 
core radius which is the region that dominate the flux. Relative to 
the pion decay emission, the sIC is suppressed by a factor of ~ 60 
which can be understood by considering hadronic decay physics 
and the fact that the CR energy spectrum, E 2 f(E), is decreasing as 
a function of proton energy (see equation[T9t. Even though the pIC 
emission component is sub-dominant, it shows a rather flat spec- 
tral index. This implies only a few strong shocks that are respon- 
sible for the electron acceleration. These merging shock waves are 
traversing the cooling core region in the cluster center. We caution 
the reader that the over-cooling of the cluster centers in our simula- 
tions possibly overestimates the true shock strengths and numbers 
which also results in an artificially enhanced pIC emission. At high 
energies, the electron cooling time is smaller than the time scale 
for diffusive shock acceleration which causes an exponential cutoff 
in the electron spectrum which is passed on to the pIC spectrum. 
The energy scale of the cutoff £ic, C ui (combining equations [6] and 
1121 ) scales with the magnetic field which causes the low magnetic 
field model to turn down faster than the large magnetic field model. 
Note that the second cutoff in the figure for the small central mag- 
netic fields is caused by a small fraction of the particles that have 
unusually high electron energy cutoff. The lower panel shows the 
photon spectral index defined in equation i I St . where T ^ 1.3 for 
the 7r°-decay emission after the pion bump that slowly flattens out 
with energy. For both the pIC and sIC emission the photon spectral 
index T ^ 1.1 above the MeV regime up to at about 100 GeV. The 
reason for the flat spectra is that the pIC emission is generated by 
a few strong shocks in the center, while the sIC emission is caused 
by protons in the flat high energy part of CR spectrum. 

Now we turn to the y-ray spectra in the WHIM which are 
shown in the right panel of Fig. [5] for our g72a cluster. Here we 
define the WHIM by the emission in the region < R < 3 R V1I 
as seen in 2D projection of the cluster. We see that the pion decay 
and sIC are suppressed by a factor that is larger than 10 compared 
to the flux within the core region since the emission correlates with 
the CR- and gas densities. The suppression of y-ray flux emitted by 
the intergalactic medium is expected to be much greater due to the 
large density decrease. The presence of small groups in the super- 
cluster region with densities that are much larger than the average 
density in the WHIM partially counteracts the flux suppression. 
Contrary to the pion decay and sIC component, the pIC emission is 
boosted by a factor of a few compared to the center because of the 
larger spatial region in the WHIM that contains a greater number of 
shocks. This leads to comparable flux from the pIC and pion decay 
emission above the energy of the pion bump in the super-cluster re- 
gion. Note however, that this flux is still sub-dominant compared to 
the pion decay flux emitted by the cluster core region. The different 
central magnetic fields do not play any significant role in the power- 
law regime in the WHIM since B <k Bcmb, which implies that the 
CRes mainly cool through IC emission. However, note that the pIC 
cutoff is shifted towards lower energies for these smaller magnetic 
fields since £ ICiBIH oc B/(B 2 + S^ MB ) oc B/Bj; MB (as derived by 
combining equations l6land lBl31 >. In the lower panel we show that 
the photon spectral index is steeper in the WHIM for all three emis- 
sion components compared to the core region. The photon index is 
about 1.3 for both pIC and sIC below 10 GeV and about 1.4 for the 
pion decay above the energy of the pion bump. The reason for the 
steeper pIC is because most of the energy that is injected into pri- 
mary electrons comes from multiple intermediate-strength shocks 
(accretion shocks), while in the cluster center the pIC emission is 
build up from a few strong shocks in the over-cooled center (merger 
shocks). The steeper sIC and pion decay photon indices are caused 
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Figure 6. Spectral universality across our cluster sample: we show the me- 
dian CR proton spectrum (solid line) as a function of dimensionless CR 
momentum p = P p /m p c for our sample of 14 galaxy clusters in Table [TJ 
together with the 68 percentiles (dotted lines). The CR spectrum of every 
cluster, f v (p), has been obtained by volume weighting the individual spec- 
tra of our SPH particles and normalizing them at p = 1. Note the small 
variance between different clusters below the diffusive break, which indi- 
cates a universal CR spectrum for galaxy clusters. The large scatter in the 
p = 10 6 - 10 8 momentum regime is a consequence of the strong radial 
dependence of the diffusive break. The lower panel shows the CR spectral 
index a = -d log //(d log p) of the cluster sample. 

by the slightly steeper CR spectrum present in the WHIM of the 
g72a cluster (cf. Fig.|7J. 



4 THE CR PROTON DISTRIBUTION 

In this section we investigate the CR proton spectrum that we ob- 
tain from our simulations and discuss the relevant physics that gives 
rise to it. We explore the variance of the spectrum across our cluster 
sample and within individual clusters and show that it obeys a uni- 
versal spectral shape. In addition, we study the spatial profile of the 
CRs within a cluster as well as across our cluster sample and find 
it to be approximately universal. This universal behavior enables 
us to construct a semi-analytic CR spectrum and to compute the 
y-ray spectrum as well as other secondary decay spectra (electrons, 
neutrinos) semi-analytically. 

4.1 A universal CR spectrum 

In Fig.[6]we show the median CR spectrum of all 14 galaxy clusters 
as well as the associated spectral index alpha. The CR spectrum 
of every cluster, / v (p), has been obtained by volume weighting the 
individual spectra of our SPH particles which have been normalized 
at the dimensionless proton momentum p = P p /m p c = 1. Before 
discussing the spectral shape we note that it shows a remarkably 
small variance across our cluster sample which indicates a universal 
CR spectrum for galaxy clusters. There are three important features 
in the spectra. 

(i) The spectrum shows a low-momentum cutoff due to ef- 
ficient Coulomb cooling at these low momenta with p < 1: 
the CR energy is transferred into the thermal energy reservoir 
through individual electron scatterings in the Coulomb field of the 
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CR particle (Gould 1 19721) . The Coulomb time scale of a mono- 
energetic CR population is very short, t Cou i = £cr/^cr.Cou1 - 

0. 03 Gyr x (p/0. 1) 3 (« e /10~ 3 cm~ 3 ) , where we show a momentum 
scaling that is valid only for the relevant non-relativistic regime. 
The Coulomb time scale for a power-law population of CRs can 
be significantly longer, t Coii1 = s C r(C, q, ar)/s C R(C, q, a) Cou i - 

1 Gyr x (ra e /10~ 3 cnT 3 )~', where we assumed a low -momentum 
cutoff q = 0.1 and a = 2.5 dEnfllin et al.1 120071) . This, how- 
ever, is still short compared to the dynamical time scale r dyn ~ 

2 Gyr x (« e /10~ 3 cnr 3 ) -1 ' 2 , Hence, we expect the formation of an 
equilibrium cutoff of our CR spectrum around q = 0. 1 ... 1 for 
the typical number densities n c ~ (10~ 3 . . . 10~ 2 ) cm -3 that we en- 
counter at the cluster cores. Note that in the presence of cooling 
processes only, the equilibrium cutoff q is determined from the 
competition between Coulomb cooling and hadronic losses and 
con verges around q = 1.685 if the cooling time is sufficiently short 
(see ljubelgas et al.l d2008h for a detailed discussion). The reason for 
that is that Coulomb cooling shifts the cutoff to higher momenta, as 
the CRs with low momenta are transferred to the thermal reser- 
voir. At high momenta, the cooling time due to hadronic interac- 
tions is shorter than the Coulomb cooling time. Hadronic cooling 
effectively removes the CRs with high energy and moves the cutoff 
towards lower momenta. 

(ii) In the momentum range between p = 1 - 10 6 , the spec- 
trum has a concave shape in double-logarithmic representation, 

1. e. it is a decreasing function with energy in the GeV/TeV en- 
ergy regime. This is quantified by the momentum dependent spec- 
tral index (shown in the lower panel of Fig. [6j( which ranges from 
a ~ 2.5 at energies above a GeV to a ~ 2.2 around 100 TeV. This 
spectral shape is a consequence of the cosmologic al Mach number 
distri bution that is mapped onto the CR spectrum jPfrommer et al.l 
2006). This mapping depends on the shock acceleration efficiency 
as a function of shock strength as well as on the property of the 
transport of CRs into galaxy clusters. Nevertheless, we can easily 
understand the qualitative features: the typical shocks responsible 
for CR acceleration are stronger at higher redshift since they en- 
counter the cold unshocked inter-galactic medium|3 This implies 
the build-up of a hard CR population. Since the forming objects 
have been smaller in a hierarchically growing Universe, their grav- 
ity sources smaller accretion velocities which results in smaller 
shock velocities, u. Hence the energy flux through the shock sur- 
faces, Efas/R 2 oc pv 3 , that will be dissipated is much smaller than 
for shocks today. This causes a lower normalization of this hard 
CR population. With increasing cosmic time, more energy is dissi- 
pated in weaker shocks which results in a softer injection spectrum. 
Despite the lower acceleration efficiency, the normalization of the 
injected CR population is larger that that of the older flat CR pop- 
ulation which yields an overall concave spectral curvature. We will 
study the details of the CR acceleration and transport that leads to 



7 Note that after re-ionization, the UV background sets a temperature floor 
by photo-ionizing the intergalactic medium (IGM) to temperatures of about 
T ~ 10 4 K except for void regions where the densities are too small so 
that the photo-io nization heating rate is smaller than the expansion rate 
iKatz et alj[l99r3) . As time progresses, formation shocks dissipate energy 
and heat the IGM to temperatures of T ~ (10 5 - 10 7 ) K which is even 
heated further as it is accreted onto clusters. Since the sound-speed is pro- 
portional to r 5 , we expect a higher sound-speed in (pre-)collapsed objects 
at low redshifts which results in weaker shocks at later times. Th is is quan- 
tified in the redshift evolu tion of the Mach number distribution iRvu et al] 
l2003l ; |Pfrommer et alj20"53) that also confirms the argument given above. 
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Figure 7. Spectral universality within a cluster - variance and radial depen- 
dence: we show the CR proton spectrum of a large merging cluster, g72a, 
as a function of dimensionless CR momentum p = P p /m p c. Top: median 
of the volume weighted and normalized CR spectrum f v (p)p 2 (solid) and 
particle-by-particle variance as indicated by the 68 percentiles (dotted line). 
Bottom: volume weighted CR proton spectrum (median and 68 percentiles) 
for different radial bins. Those are represented by different colors: core re- 
gion of the intra-cluster medium in red, the intermediate cluster scales in 
orange, the periphery of the ICM in green, and the WHIM in blue. The 
spectrum for the different radial bins has been offset vertically by a constant 
factor of four, with decreasing amplitude for increasing radius. The shift of 
the cutoff ijtoa higher momentum for smaller radius is due to the smaller 
CR cooling time scales in denser regions. Note that the cosmic ray spectrum 
inside the virial radius is almost independent of radius. This behavior is not 
only observed for this cluster, but persistent across our cluster sample. 

that particular spectrum in a forthcoming paper (Pinzke & Pfrom- 
mer, in prep.). 

(iii) There is a diffusive break in the spectrum at high momenta 
where the spectral index steepens by 0.3. The CRs above these en- 
ergies can diffusively escape from the cluster within a Hubble time. 
The particular value of the steepening assumes that the CRs scat- 
ter off Kolmogorov turbulence. Using twice the virial radius of each 
cluster, we find that the diffusive break varies between the momenta 
p = 10 6 - 10 s , dependent on the characteristic size of a cluster 
(equation [3j. 

We now turn to the question on how universal is our CR spec- 
trum within a cluster. In Fig. [7] we find that the intrinsic scatter of 
our CR spectra within a cluster is larger than the scatter among 
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clusters. The reason being that formation shocks are intermittent 
as mass is accreted in clumps an not continuously. Hence there is 
a high intrinsic variance of the CR spectrum for similar fluid ele- 
ments that end up in the same galaxy cluster. However, averaging 
over all fluid elements that accrete onto a galaxy cluster results in a 
very similar spectrum since different galaxy cluster experience on 
average the same formation history^ We note that the spectral vari- 
ance of g72a is representative for all the CR spectra in our sample. 

We study the radial dependence of the CR spectrum for g72a 
at the bottom of Fig. [7] Inside the cluster, the spectral shape does 
not strongly depend on the radius. This is a crucial finding as it 
enables us to separate the spectral and the spatial part of the CR 
distribution. The level of particle-by-particle variance is similar to 
that of the total cluster spectrum. Also noticeable is the increasing 
low-momentum cutoff q as we approach the denser cluster center. 
This is due to enhanced Coulomb losses and to a lesser extent in- 
creased adiabatic compression that the CR distribution experienced 
when it was transported there. Outside the cluster g72a, for radii 
R > R vn , we observe a considerably steeper CR spectrum at CR 
energies of E i 1 TeV compared to the cluster region. We note that 
this behavior is not universal and we observe a large scatter of the 
CR spectral indices among our cluster sample at these large radii. 
This behavior might be caused by the recent dynamical activity of 
the cluster under consideration but a detailed characterization goes 
beyond the scope of this work and will be postponed (Pinzke & 
Pfrommer, in prep.). 

4.2 The spatial distribution of CRs 

We have shown that the CR spectrum is separable in a spectral and 
spatial part. To this end, we quantify the spatial part of the CR spec- 
trum by taking the volume weighted CR spectrum in each radial bin 
i (see Section lATl in the appendix for derivation), 



f v (Ri) = </> v (p = 1, J? f ) = C y (Ri) = C M (R,) ' 



(20) 



Here we assume q < 1, the subscripts M and V denote mass- and 
volume weighted quantities, respectively, and we introduced the 
dimensionless normalization of the CR spectrum, 



C(R) = rru 



C(R) 



m p (a- 1) q a 



, n CR (R) 



(21) 



p p(R) " pv "- 1 p(R) 

where iiqr is the CR number density. We now have to take into 
account physical effects that shape the profile of C. Those include 
acceleration of CRs at structure formation shocks, the subsequent 
adiabatic transport of CRs during the formation of the halos as well 
as non-adiabatic CR cooling processes; primarily hadronic interac- 
tions. At the same time we have to consider the assembly of the 
thermal plasma and CRs within the framework of structure forma- 
tion that is dominated by CDM. Hierarchical structure formation 
predicts a difference in the halo formation time depending on the 
halo mass, i.e. smaller halos form earlier when the average mass 
density of the Universe was higher. This leads to more concentrated 
density profiles for smaller halo masses; an eff ect that breaks th e 
scale invariance of the dark matter halo profile fehao et all 2009), 
The central core part is assembled in a regime of fast accretion 



The reason for our increased scatter at higher energies is that we nor- 
malize the spectrum of each SPH particle at p = 1 where the softest CR 
populations dominate. Thus, at high energies the variance is mainly driven 
by the scatter in the hard CR populations, which have been accelerated in 
strong shocks at higher redshifts. 



(Lu et "aT] |2006l) . This violent formation epoch should have caused 
the CRs to be adiabatically compressed. In the further evolution, 
some cluster have been able to develop a cool core which addi- 
tionally could have caused the CRs to be adiabatically contracted. 
On larger scales, the gas distribution follows that of dark matter (at 
least in the absence of violent merging events that could separate 
both components for time scales that are of order of the dynami- 
cal time). On these scales, the CR numbe r density roughly traces 
the gas distribution ( Pfro mmer et"ail 2008). Overall, we expect the 
spatial CR density profile relative to that of the gas density to scale 
with the cluster mass. If non-adiabatic CR transport processes have 
a sufficiently weak impact, these considerations would predict flat- 
ter C-profiles in larger halos as these halos are less concentrated. 

Our goal is to characterize the trend of C-profiles with cluster 
mass and to test whether we can dissect a universal spatial CR pro- 
file. To this end, we adopt a phenomenological profile that allows 
for enough freedom to capture these features as accurately as pos- 
sible. Hence, we parametrize C(R) with shape parameters that in- 
clude a flat central region given by C cmta , a transition region where 
the location is denoted by R lmns and the steepness by /?, and finally 
a flat outer region denoted by C VK , through the equation 



Cm(R) - (C v i r - Ccenter) 1 + 



^tran 



-0 



(22) 



The core regions in our radiative simulations show too much 
cooling, and we possibly lack of important CR physics such as 
anisotropic CR diffusion that could smooth out any strong inho- 
mogeneity at the center. Hence we adopt a conservative limit of the 
central profile value of C centcr = 5 x 10~ 7 . 

In the top panel of Fig. [8] we show the mean of the Cm pro- 
files from the cluster simulation sample. We subdivide or sample 
in two different mass intervals: large mass clusters (top red), and 
small mass clusters (bottom blue), in the mass range 1 x 10 15 < 
M vir /M < 3 x 10 15 , and 7 x 10 13 < M vir /M < 4 x 10 14 , re- 
spectively^] The error-bars show the standard deviation from sam- 
ple mean and the solid lines the best fit that will be discussed in 
more detail in Section 1431 The Cm profile of our large mass clus- 
ters is almost flat and shows only a very weak radial dependence. 
In contrast, the C profile of our small clusters has a rather steep 
and long transition region and is increasing towards the center. The 
difference in normalization, transition width, and transition radius 
between low mass and large mass clusters indeed suggests that Cm 
scales with the cluster mass in a way that we anticipated. We quan- 
tify the mass scaling of these shape parameters by a power-law fit to 
the small and large clusters in the three lower panels of Fig. [8] the 
value of Cm in the asymptotically flat regime in the cluster periph- 
ery, C v i r (top); the transition radius, /? tnms (middle); and the steep- 
ness of the transition, [} (bottom). As expected, we find that all three 
quantities increase slowly with mass. We additionally show the val- 
ues of Cm at R vil for each cluster (black crosses). Within the scatter, 
these values are consistent with the mass scaling found by the pro- 
file fits in our two cluster mass bins. We have shown that there 
exists an almost universal spatial CR profile after taking into ac- 
count the weak trends of the C profile with cluster mass. Note that 
the particle-by-particle variance of the spatial CR profile within a 
cluster (that we address in Section |A3l in the appendix) additionally 
supports the conclusion of a universal spatial CR profile. 



9 Note that we exclude the intermediate mass cluster gld, since it had a 
very recent merger which resulted in a double-cored system with a non- 



spherical CR profile. 
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Figure 8. The top panel shows the profile of the dimensionless normal- 
ization of the CR spectrum, Cm . We show the mean Cm and the standard 
deviation across our cluster sample which has been subdivided into two dif- 
ferent mass intervals: large- (red), and low-mass clusters (blue) representing 
the mass range 1 X 10 15 < M vil -/M Q < 3 X 10 15 , and 7 X 10 13 < M vil -/M < 
4 X 10 14 . The solid lines show the best fit to equation 1221 . The lower three 
panels show the mass dependence of the quantities which parametrize Cm 
for low mass clusters (blue circles) and large mass clusters (red circles). 
The top small panel shows the asymptotic Cm for large radii (C v i r ), where 
each cross shows Cm at i? v j, for each cluster. The middle panel shows the 
transition radius R t rans. and the bottom panel shows the inverse transition 
width denoted by fi. 



4.3 A semi-analytic model for the spatial and spectral CR 
distribution 

In our simulations we use a CR spectral description which follows 
five different CR proton populations; each being represented by a 
single power-law. The CR populations are chosen such that we ac- 
curately can capture the total CR spectrum in the entire momentum 
space (a convergence study on the number of CR populations is 
presented in the appendix lA2t . We want to model this spectrum an- 
alytically with as few CR populations as possible, but at the same 
time preserve the functional form of a power-law of each popula- 
tion. In that way we can easily obtain the total CR spectrum by 
superposition and apply a simple analytic formula to estimate the 
induced radiative processes. 

In detail, we use a total CR proton spectrum f v (p,R) that is 
obtained by summing over the individual CR populations fi (equa- 
tion^; each being a power-law in particle momentum with the total 




Figure 9. Cosmic ray proton spectrum as a function of dimensionless CR 
momentum p = P p /m p c for our sample of 14 galaxy clusters (TablefT}. The 
main panel shows median of the cosmic ray spectrum f y (p)p 2 across our 
cluster sample. The spectrum of each cluster has been normalized at p = 1 
(black crosses). The blue solid line shows the best fit triple power-law to 
the simulation data. The light blue area shows the 68 percentile of the data 
points across our cluster sample. The bottom panel shows the difference be- 
tween the relative fit and the simulation data (blue solid line) which amounts 
to less than 2 percent. 



CR amplitude C V (R), 

fy(P, R) = C v (i?) g(£p,max) Dp(P, Pbreak, ?) ^ \ P~ 



(23) 



where we assume universal spectral shape throughout the cluster. 
The normalization of / v depends on maximum shock acceleration 
efficiency f pjnax , where g(0.5) = 1 and g(£ p , max < 0.5) < 1 (func- 
tional will be studied in Pfrommer in prep.). Denoting the ampli- 
tude of each CR population by Cj(R), where the number of spectral 
bins might be different from the five chosen in the simulations, we 
construct the relative normalization for each CR population 



Cj(R) 



(24) 



Here we have assumed that Cj(R)/C(R) is only a weak function of 
radius. This is explicitly shown in Fig. [7] where the functional form 
is almost independent of the radius. The two energy breaks in the 
CR spectrum are represented by 



D p (p,p bxak ,q) = 



1 + , 



\ Pbreak / 



(25) 



The first term in equation {25} ensures the low-momentum slope 
~ p 2 as appropriate from th e phase space volume that is populated 
by CRs jEnBlin et alj|2007l) and the last term accounts for diffu- 
sive losses of CRs that steepen the spectrum by 5 = 1/3 assum- 
ing Kolmogorov turbulence. The shape parameter fi determines the 
size of the transition in momentum space. Our choice of fi = 3 
ensures a fast transition from one regime to the other. The low- 
momentum break q is determined to q = 0.3 from Fig. [9] and the 
high-momentum break pbreak is derived from equation l|3j- 

To capture the spectral universality of our cluster sample, we 
fit in Fig. [9] the median of our cluster sample of the CR spec- 



© 2008 RAS, MNRAS 000,[T]f33] 



Simulating the gamma-ray emission from galaxy clusters 15 



trum with a triple power-law function. Because we normalize the 
spectrum at p = 1, this ensures that the normalized spectrum 
fv(p,R)/f v (l,R) = £; A,p~ a ' becomes independent of radius (cf. 
equation [23}, The triple power-law fit represented by the blue line 
in the upper panel of Fig.[9]resulted in 



cluster simulations 



A = (0.767, 0. 143, 0.0975) , and a = (2.55,2.3,2.15) . 



(26) 



The data from the simulation is denoted by black crosses, together 
with the 68 percentiles spread shown by the light blue area. In the 
lower panel, we show the relative difference between the simulation 
and the fit which indicates a difference of less than two percent from 
the GeV to PeV energy regime. 

The spatial part of the CR spectrum is derived from the fit to 
the mean C M in the top panel of Fig.(8] The solid lines show the best 
fit to Cm using equation &221 . We find that the central value C cen ter = 
5 x 10~ 7 is the most conservative value that still provides a good fit 
for both mass intervals. Note that there is a larg e uncertainty in 
this value due to insufficient data in the centeo, implying that it 
should be treated as a lower limit. However, the gamma-ray flux 
depends only weakly on the exact value of C cen ter since most of 
the flux reside from the region outside the transition region. The 
mass dependence of Cm is quantified in the three lower panels of 
Fig. [8] where we fit a power-law in mass for the normalization C v j, 
(top), the transition radius ^ tlans (middle), and the steepness of the 
transition region /? (bottom). We obtain the following relations, 

C vir = 1.7 x 10~ 7 x(M vil ./10 I5 M o ) - 51 , (27) 
R tmm = 0.021 J? vir x(M vil ./10 15 M G ) - 39 , (28) 
B = 1.04 x(M vir /lO 15 M ) 015 . (29) 
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Figure 10. Schematic overview of how our semi-analytic framework re- 
lates to the simulated CR and y-ray quantities. From the cluster simu- 
lations we derive the CR spectrum for each SPH particle, f(p,R). Go- 
ing clockwise, the integrated y-ray production rat e for each SPH particl e 
A n o_ ysim (E,R) is derived within the framework of IPfrommer et all ( 2008). 



The radial profile of the integrated y-ray production rate, A^ 



-CE.J). 



is obtained by volume weighting this quantity of individual SPH particles. 
Returning to f(p, R) and going counter-clockwise instead, we can directly 
perform a volume-weighting of the CR spectrum in each radial bin, yield- 
ing f s i m (p,R). Fitting both the spatial and spectral part provides us with a 
semi-analytic model for the spectrum, / ma . The semi-analytic model for the 
in tegrated y-ray pr oduction rate, ^ T o_ yana , is derived within the framework 
of IPfrommer et alj 120081) . The explicit form of ^o_ yaila can be found in 
equations d30M34t . while Aq and &(E,aj) are implicitly given by equa- 
tion i32) . Comparing the flux from the simulations and our semi-analytic 
model shows good agreement for our cluster sample. 



5 SEMI-ANALYTIC MODEL FOR THE y-RAY 
EMISSION 

In this section we derive a semi-analytic formula for the integrated 
y-ray source function that is based on our semi-analytic CR dis- 
tribution. Using the gas density profile of the cluster along with 
its virial radius and mass, this formula enables us to predict the 
dominating ;r -decay induced y-ray emission from that cluster. The 
density profiles can either be inferred from simulations or from X- 
ray data. To test the semi-analytic source function, we use density 
profiles from our cluster simulations. We also make y-ray flux pre- 
dictions for the Coma and Perseus cluster using their true density 
profiles as inferred from X-ray observations. 

5.1 Schematic overview and semi-analytic formulae 

In Fig.[l0]we show a schematic overview of the simulated CR spec- 
trum and the integrated source function together with the mapping 
to our semi-analytic model. From the cluster simulations we derive 
the CR spectrum f(p, R) describing the phase space density of CRs. 
When taking the spherical (volume weighted) average of f(p, R), 
we obtain f Y (p,R). We fit the spectral and spatial part of this func- 
tion separately. The semi-analytic model for the CR distribution in 
clusters that we provide in equation ([23} can now be used to pre- 
dict the secondary radio synchrotron emission and the hadronically 

10 We cut the Cm profiles at R =s 0.03R v j r due to the bias with the over- 
cooled center. The increased scatter in the center for the small clusters is 
caused by the low statistics of the fewer number of clusters that can con- 
tribute on these small scales. With both of these considerations in mind, we 
conclude that there is a large scatter in C ccn ter- 



induced neutrino an d y-ray emission from g alaxy clusters. Follow- 
ing the formalism in lPfrommer et alj {2008), we obtain the volume 
weighted and energy integrated and omnidirectional (i.e integrated 
over the 47T solid angle) y-ray source function due to pion decaNp^l 
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1 + 



IF ■ 



(32) 



where the sum extends over our three CR populations, CmW is 
given by equations l !22t , Ml\ - {29}, while A and a are provided by 
equation J26t and we introduced an auxiliary variable p for dimen- 
sional reasons to ensure that A A is dimensionless. In equation d32b 
we have also have introduced the abbreviation 



[SAa,b)lS =% X2 (a,b)-B xl (a,b) 



(33) 



where B x (a, b) denotes the incomplete Beta-function, we have used 
that the number density of target nucleons is the sum of hydrogen 
«h and helium «He number densities, «n = «h + 4«He = p/m p . 
The shape parameter (5, = 0. 14 or 1 - 6 + 0.44 allows us to accurately 

" From here on we always imply the volume weighted and energy inte- 
grated source function when talking about A if nothing else is stated. 
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predict the emission close to the pion bump in combination with 
the effective inelastic cross-section for proton-proton interactions, 
<x p p_, =a 32(0.96 + e 442 - 2 - 4ff i) mbarn. In addition we have a term, 
similar to equation {25} that describes the diffusive CR losses due 
to escaping protons from the cluster at the equivalent photon energy 
for the break, £ y brcak . It is derived from E„ brcak of equation {3}, 

-,-1/(3/3) 



DJE y , E 



y, break 



l + 



-<y, break 



(34) 



where /? = 3. Finally we note that the semi-analytic y-ray model is 
derived within i? v jj. Applying the model to the region outside R ya , 
but within "iR yn , would increase the y-ray flux by less than 10% for 
small mass clusters and less than 30% for large mass clusters (cf. 
Table©. 

For convenience and completeness, we provide here the dif- 
fer ential y-ray source functi on for pion decay (for further details 
see lPfrommer & Enfilirj2004h . 



s n o_ y (R, E) 
s n o- y (E) 



<?(^Tp.max) ^y(^y; $y, break) 
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3iripC 



Po 



2m„ 



Ai 



2£ v 



2E V 



(35) 



(36) 



where A S (R) = A A (R) (see equation |31t . The exact shape of 
D 7 {E y , E 7t break) depends on the detailed physics of CR diffusion 
and the characteristics of turbulence and is subject to future stud- 
ies. We note, however, that the break does not interfere with the 
energy range of current y-ray observatories. 



5.2 Comparing our semi-analytic model to simulations 

To test our semi-analytic y-ray model we contrast it with numer- 
ically calculated radial profiles and spectra. In the upper panel of 
Fig. [TT] we compare the radial profile of our semi-analytic y-ray 
source function (equation[30j at 100 MeV to a numerical emission 
profile for two representative clusters, a large post-merging cluster 
g72a, and a small CC cluster g914. The numerical profile has been 
obtained by means of calculating the y-ray source function A.^)_y of 
every SPH particle and volume weighting the resulting radial pro- 
file. The overall shape of the radial profiles of A !T o_ y for different 
clusters are quite similar. This is because the main spatial depen- 
dence originates from the gas density, that enters with a square in 
A.„o_y and has a small scatter across the cluster sample. The differ- 
ent behavior in the cluster centers stems from the steeper profile 
of Cm(R) for low mass clusters. We show the difference between 
the integrated source functions from our simulations and the semi- 
analytic model in more detail in the lower panel of Fig.QT] In both 
clusters, we see an excellent agreement with differences amount- 
ing to less than 20-30 percent at any radius. These differences are 
representative for our cluster sample, which is indicated by the me- 
dian difference across our cluster sample together with the 68 per- 
centiles that show a spread similar to the difference in these two 
clusters. Since these are fluctuating differences, they average partly 
away when we perform the volume integral. Hence we obtain an 
agreement of the total flux within the virial radius - obtained di- 
rectly from the simulations and in our semi-analytic model - that 
is better than 5 per cent for the two representative clusters that are 
shown. This indicates that our semi-analytic description manages to 
capture the CR physics, that is important for y-ray emission from 
clusters, surprisingly well. 
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Figure 11. Comparison of the spatial profile of the y-ray source function 
in our simulations and our semi-analytic framework. The main panel shows 
the profile of the y-ray source function A^s_ y that results from pion decay 
emission. We compare /l^o_ y at 100 MeV for two clusters, the large merging 
g72a (red) cluster and the the small CC g914 (blue) cluster, both as a func- 
tion of R/R vn . The solid lines show the simulated integrated source func- 
tion and the dotted lines the semi-analytic model. The lower panel shows 
the difference between the semi-analytic and the simulated source func- 
tions normalized by the simulated source function. The orange lines show 
the median (solid) difference between the semi-analytic and the simulated 
source functions across our cluster sample together with the 68 percentiles 
(dotted). 

We show the y-ray spectrum of the 7r°-decay emission 
weighted by energy for g72a and g914 in the upper panel of Fig. [12] 
Both clusters have very similar spectra with the exception of the 
diffusive steepening that is inherited from the proton spectrum. 
Since the break of the proton spectrum scales as £ p _ break x Mj ir 
with the virial mass of the cluster (equation [3}, the break in the 
pion decay spectrum is reduced by a factor of 10 3 for the smaller 
cluster g914. The solid and dotted lines contrast the simulated spec- 
trum to that obtained from our semi-analytic model. The difference 
between these two approaches amounts to less than 20 percent for 
both clusters in the GeV and TeV band (shown in the lower panel). 
The flux differences between our semi-analytic model and the sim- 
ulations for individual clusters are about a factor two smaller com- 
pared to the scatter in the the mass-luminosity scaling relations for 
a given cluster (see e.g. Fig.|15ll. The reason for the more accurate 
predictions within our semi-analytic formalism is a direct conse- 
quence of the essential additional spatial information of the gas and 
CR density that we account for. 



5.3 Predicting the y-ray emission from Perseus and Coma 

Here we demonstrate how our semi-analytic formalism can be ap- 
plied to predict the y-ray flux and surface brightness from real clus- 
ters using their electron density profile as inferred from X-ray mea- 
surements. The predicted flux and surface brightness are then com- 
pared to current upper limits and previous work. 

The two clusters that we investigate are two of the 
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Table 3. Properties of the Coma and Perseus galaxy cluster. 



cluster 


z m 


D\ [) 

Jum 


vir 


M (1) 
vir 


L (1) 
XO.l-2.4 


n m 
e 






[Mpc] 


[Mpc] 


[10 14 M G ] 


[10 44 ergs"'] 


[10~ 3 electrons cm" 3 ] 


Coma 


0.0232 


101 


2.3 


13.8 


4.0 


r , n -l.l25 
3.4 X [1 + (i?/294kpc) 2 ] 


Perseus 


0.0183 


77.7 


1.9 


7.71 


8.3 


46x [l +(i?/57kpc) 2 ] _LS +4.79X [l + (R/200kpc) 2 ]~°' 87 



Notes: 

(1) Data for redshift (z), luminosity distance (£*i um ), /? v i r , M yK , and X-ray 0.1-2.4 keV luminosity (Lx,o. 1-2.4) are taken from Reiprich & Bohringer (2002) and 
rescale d to our assumed cosmolog y (Section f2.U . The electron number density, n e , for Coma and Perseus are inferred from X-ray observations bv lBriel et alj 
1 19921) and lChurazov et al l 120031) . respectively. 

Table 4. y-ray flux comparison of the Coma and Perseus galaxy cluster. 







E 7 
[GeV] 


Experiment 


9 

[deg] 


5>,UL(> Ey) 

[phcirT 2 s~'] 


[phcm -2 s" 1 ] 




GeV-regime 


Coma 


100 MeV 


EGRET* 1 ' 


5.8< 4 ' 


3.81 X 10~ 8 


4.20 x 10~ 9 


9.1 




Perseus 


100 MeV 


egret' 1 ' 


5.8< 4 ' 


3.7 X 10~ 8 


1.46 x 10~ 8 


2.5 


TeV-regime 


Coma 


1 TeV 


HESS (2) 


0.2 


1.1 xlO~ 12 


4.86 x 10~ 14 


23.0 




Perseus 


1 TeV 


Magic' 3 ' 


0.15 


4.7 x 10~ 13 


1.75 x 10~ 13 


2.7 




Perseus 


100 GeV 


Magic' 3 ' 


0.15 


6.55 x 10~ 12 


3.2 x 10~ 12 


2.0 



N otes: 

m lReimer et al.1 120031) . (2) lAharonian et alj 120091) . ffl lAleksic et alj feoiot) , (4) Since the flux is dominated by the region inside R vir we use R e = R vil . 
(5) y-ray fluxes are obtained within our semi-analytic model that is based on our conservative model. The y-ray flux level depends on the maximum shock 
acceleration efficiency of CRs, f p . raax , for which we assume an optimistic value of 0.5. Smaller efficiencies imply smaller fluxes. 

5.3.1 Comparison to upper limits in the GeV/TeV regime 



We calculate the y-ray fluxes for Coma and Perseus above 100 MeV 
and 100 GeV that is emitted within R e = _R v j, (see Table|4]l. Com- 
paring the 100 MeV-flux in our model to the EGRET upper limit on 
the Coma (Perseus) cluster (Rei mer et al.l2003T) . we find that it falls 
short of our semi-analytic prediction by a factor of about 9 (2.5). 
With the two-year data by Fermi, this upper limit on Coma will 
improve considerably and is expected to become competitive with 
our predictions. With this at hand, we will be able to put impor- 
tant constraints on the adopted CR physics in our simulations. In 
particular, we can test our assumptions about the maximum shock 
acceleration efficiency at structure formation shocks. Since Fermi 
detected y-rays from the central cD galaxy in Perseus, NGC1275, 
at a level that is about five times higher than the EGRET upper lim- 
its, this indi cates that the sou rce is variable on time scales of years 
to decades jAbdo et alj[2009h . Hence it restricts and complicates 
the detectability of the extended pion-day emission that might be 
buried underneath. 

In the TeV regime, we integrate our model prediction within 
a solid angle that is comparable to the point-sp read function of 
IACT s. The best current upper limit for Coma jAharonian et al.l 
2009) falls short of our semi-analytic prediction by a factor 20. 
The best curr ent upper flux lim it for Perseus (using a spectral 
index of -2.2, Aleksi c et a l. 2010) is only a factor of two larger 
than our flux prediction and clearly within reach of future deeper 
TeV observations. This demonstrates the huge potential of nearby 
CC clusters to detect non- thermal y-ray emission as suggested by 
Pfrommer & Enfilinl l [2004l) . 



brightest X-ray clusters in the extended HIFLUGCS catalogue 
(Reiprich & Bohringer! l2002h - a sample of the brightest X-ray 
clusters observed by ROSAT - namely Coma and Perseus. Both 
Coma and Perseus are well studied clusters, where Coma is a large 
post-merging cluster while Perseus is a somewhat smaller clus- 
ter that hosts a massive cooling flow and is the brightest X-ray 
cluster known Edge et alj 1 1992). In Table [3] we show the data 
for respective cluster. Using the electron number density profile, 
we can calculate the gas density profile of the cluster through 
p = m p n s /(XuX c ). Here denote Xn = 0.76 the primordial hy- 
drogen (H) mass fraction, and the ratio of electron and hydrogen 
number densities in the fully ionized ICM is given by X e = 1.157. 
For sufficiently small angular scales, the y-ray flux within the ra- 
dius R s = 9D\ um of a disk of angular radius 8 (measured in radians 
and centered at R = 0) is calculated by 



1 r R « 

T^ y (>E y ) = — 2ndR ± R ± S^_ y (R ± ,E y ), (37) 

D ium JO 

S ^ {R± , Ey) = 2^ r«^4^, (38) 



An 



Po V* 2 - Rl ' 



where A R o_ y (E y ) is given by equation J32b and we introduced 
the definition for the y-ray surface brightness S R o_ y in the last 
step. As a consistency check, we compare the flux ratios from X- 
rays 7 r x-iay.Pei.scu S /9 r x-ra y ,coma - 3.51 to that of y-ray s as predicted 
by our semi-analytic model, 7v>_ r ,p crsclls (> 100MeV)/yyi_ 7]Coma (> 
100 MeV) ^ 3.47 and find excellent agreement within one percent. 
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Figure 12. Comparison of the y-ray spectrum in our simulations and our 
semi-analytic framework. The main panel shows the y-ray number flux 
from pion decay weighted by the photon energy. We compare two clusters, 
the large merging g72a (red) cluster and the small CC g914 (blue) cluster, 
both as a function of photon energy. The solid lines show the simulated flux 
and the dotted lines the semi-analytic model. The lower panel shows the 
difference between the semi-analytic and the simulated flux normalized by 
the simulated flux. 

5.3.2 Surface brightness profile 

To explain the large difference in flux between Coma and Perseus 
we show the surface brightness as a function of the viewing angle 6 
in Fig.[T3] The dense cooling core region of Perseus provides ample 
target material for two-body interactions such as bremsstrahlung or 
hadronic CR interactions which boosts the luminosities likewise. 
This results in a large increase in flux compared to the average 
cluster of similar mass that are characterized by mass-luminosity 
scaling relations (see e.g. Table [5}. Assigning a flux to a cluster 
that is consistent with the scaling relations should be a very conser- 
vative approach for CC clusters. The half flux radius for Perseus 
(Coma) is f?HF = 0.11 deg (0.18 deg). It is shown with dotted 
lines in Fig.[T3] Since these f? H F values are comparable to the angu- 
lar scale of the point spread functions of IACTs (0.1-0.2 degrees) 
both clusters are suitable candidate sources. The dashed and dotted- 
dash ed lines are obtained by u sing the semi-analytic formalism 
from lKushnir & Waxmanl J2009h to predict the surface brightness 
from the Coma cluster above the energy 100 MeV and 100 GeV, 
respectively. See Section [8T| for a detailed comparison to their re- 
sult. 



5.4 CR-to-thermal pressure and temperature profile in 
Perseus and Coma 

A quantity that is of great theoretical interest is the CR pressure rel- 
ative to the thermal pressure, X C r = Pen/Ah as it directly assesses 
the CR bias of hydrostatic cluster masses where the CR pressure en- 



10 J F 




in- 7 1 , , , , i , , , 

1 10 
Angular size G [arcmin] 



Figure 13. The y-ray brightness profiles derived with our semi-analytic for- 
malism (solid line) as a function of the angular size on the sky 8. The left 
axis indicates the emission for E y > 100 MeV, and the right axis shows the 
emission for E y > 100 GeV. We compare the emission from the Perseus 
galaxy cluster (red solid line) to the Coma galaxy cluster (blue solid line) 
while employing the X-ray density profile for each cluster, respectively. 
The high surface brightness in the cool core cluster Perseus is a result of 
the high central gas densities in this cluster. The dotted line shows the ra- 
dius from within half the flux originates. It is clearly smaller for Perseus 
due to the steep central density profile and the larger transition region in 
Cm for smaller clusters. The dashed and dash-dotted lines sho w the Coma 
surface brightness as predicted from Kushnir & Waxman ( 2009) using their 
assumed parameters above the energy 100 MeV and 100 GeV, respectively. 

ters in the equation of motion. Since X C r x W th/Pih = 1/kT in the 
external cluster regions, we have to accurately model the tempera- 
ture profile of our clusters. Note, however, that this relation should 
only hold for regions with long thermal cooling times compared 
to the dynamical time scale. In particular it breaks down towards 
the center of a cooling flow cluster where the thermal gas cools 
on a shorter time scale such that the forming cooling flow causes 
adiabatic contraction of the CR population. We model the central 
regions of Coma and Perseus according to X-ray observations by 
iBriel et al.l d 19921) and lChurazov et al.l a2003h . respectively. These 
observations are not sensitive to the outer temperature profile due 
to the high particle background for XMM-Newton and Chandra!* 2 ! 
X-ray observations of somewhat more distant cluster sample show 
a universal declini ng temperature profil e outside the cooling core 
region up to R 50O jVikhlinin et alj|2005l) . We model this behavior 
of the temperature profiles towards the cluster periphery acco rding 
to cosmological cluster simulations bv lPfrommer et alld2007t) and 
obtain a function 7" cxt (A) that accounts for the decreasing temper- 
ature profile outside the core region. It is unity in the center and 
then smoothly decreases until the virial radius beyond which we 
expect the spherical approximation to break down and where the 
cluster accretion shocks should introduce breaks in the tempera- 
ture profile^ This function can be multiplied to existing central 

12 We note that Suzaku has an approximately ten times lower background 
due to its low-Earth orbit that could enable such observations in the X-rays. 
Alternatively, combining the X-ray surface brightness with future Sunyaev- 
Zel'dovich measurements of these nearby clusters should in principle allow 
the derivation of the temperature profile of the entire cluster. 

13 The mean tempe rature profile of the radiative simulations by 
IPfrommer et alj 120071) increases towards smaller radii until R ~ 0. lR v ir 
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Figure 14. Radial profiles of the gas density p, the temperature kT, and the 
CR-to-thermal pressure (Xcr> = (PcR)/(Pth)- We compare the cool core 
cluster Perseus (red) to the non-cool core cluster Coma (blue ). The density 
and and temperature pr ofiles are taken from the X-ray data iChurazov et alj 
|2003| ; iBriel et alj| 19921) while we remodel the external temperature profile 
to bring it into agreement with cosmological cluster simulations as well as 
higher redshift X-ray observations. We obtain the CR-to-thermal pressure 
profile from our semi-analytic modeling in combinations with the presented 
gas profiles. 



where R vir can be obtained from Table[3] The density and tempera- 
ture profiles for Coma and Perseus are shown in the first two panels 
of Fig.[l4] 

In the previous sections, we have seen that the y-ray surface 
brightness is a radially declining function and so is the CR pres- 
sure. In contrast, outside the central cooling core regions, the CR- 
to-thermal pressure Xqr increases with radius as can be seen in the 
bottom panel of Fig. [14] This increase is entirely driven by the de- 
creasing temperature profile. The two X C R-profiles for Coma and 
Perseus show our expectations for a typical non-CC and CC clus- 
ter. The XcR-profile in a CC cluster shows an additional enhance- 
ment towards the cluster center which results from the centrally 
enhanced CR number density due to adiabatic contraction during 
the formation of the cooling flow. During this process, the thermal 
gas cools on a short time scale compared to that of the CRs which 
causes an increase in density and hence adiabatic compression of 
the CRs. We note that the overall normalization of Xqr depends 
on the normalization of the CR distribution that itself is set by the 
maximum shock acceleration efficiency. The overall shape of X C r, 
however, should remai n invariant since CRs are adiabatically trans- 
ported into the cluster jPfrommer et alj|2007l , Pfrommer in prep.). 

While the overall characteristics of X C r in our semi-analytical 
model is similar to t hat obtained in cosmological simulations 
jPfrommer et alj|2007l) , there are some noticeable differences par- 
ticularly in the central cooling core region around the cD galaxy 
of these clusters. Again this can be traced back to known short- 
comings of modeling the physics in the central regions cor- 
rectly in current simulations such as to include AGN feed- 
back and anisotropic conduction in combination with magneto- 
hydrodynamics. This also leads to different simulated temperature 
profiles in the center compared to those inferred from X-ray obser- 
vations and explains the discrepancy in the XcR-profiles. The vol- 
ume average of the CR-to-thermal pressure for Coma and Perseus 
is (X C r) = (Pcr) I (Pth) = 0.02, dominated by the region around the 
virial radius. These values assume an optimistic saturation value of 
the shock acceleration efficiency of f max = 0.5 and decrease accord- 
ingly if this value is not realized at the relevant structure formation 
shocks responsible for the CRs in clusters. 



6 HIGH-ENERGY SCALING RELATIONS 



temperature profiles to yield the temperature profiles for Coma and 
Perseus, 

7 k eVx A +(R/7 _ lkpC)3 . r^), (39) 



^TperseusW 

kT Cma (R) 
T exl (R) 



2.3 + (tf/71kpc) 3 
8.25 keV x T cM (R), 

2-.-0.32 
/ K \ ' 

1 + 



0.2 R v 



(40) 
(41) 



due to adiabatic compression of the gas and starts to drop sharply towards 
smaller radii where radiative cooling causes the temperature to decline. 
This behavior qualitatively matches the results of l ow mass clusters from 
a Chandra sample of nearby relaxed galaxy clusters jVikhlinin et"ai]|2006t) 
whereas the temperature maximum for more massive clusters seems to shift 
to somewhat larger radii around R ~ 0.2i? v ir. Hence we adopted this larger 
value as the transition radius for Perseus and Coma. We note that the result- 
ing profiles are consistent with the o bserved central temperature profiles 
(Briel et al. 2001 : Churaz ov et alj2003h as well as the mosaiced entire tem- 
perature profile of Perseus beyond R200 (S. Allen, private communication). 



We now discuss the scaling relations of the numerical y-ray emis- 
sion from clusters and analyze their dependence on dynamical 
state, emission region and address the bias of galaxies to the total 
luminosity. The cluster scaling relations are derived by integrating 
the surface brightness map of each cluster. By fitting the total y-ray 
emission of the 14 clusters in our cluster sample with a power-law, 
we determine the mass-to-luminosity scaling. 

In the preceding sections we have shown that the pion decay 
emission dominates the total y-ray emission but we have not ad- 
dressed the question, which radii contribute most to the luminosity? 
To answer this, we have to consider the y-ray luminosity resulting 
from pion decay within a radius R, 
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Jo 



£j>JR)°c dR'R' 2 C M (R') P (R'Y 



f 

Jo 



dR' R' 2 p(R') 2 . 



(42) 



For the purpose of this simple argument, we neglected the very 
weak spatial dependence of the CR distribution which is described 
by Cm(R)- The y-ray luminosity -£^o_ r is dominated by the region 
around the scale radius R s which can be seen by considering the 
contribution to -C n o_ y per logarithmic radius, 
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Figure 15. Cluster scaling relations of the y-ray luminosity for the energy 
regimes corresponding to the Fermi y-ray space telescope and imaging air 
Cerenkov telescopes. Show is the contribution from individual y-ray emis- 
sion components within one virial radius (pIC - red squares, sIC - green 
circles, n° decay - blue triangles) to the total cluster luminosities (black 
crosses) and mass to luminosity scaling (black solid line) at the energies 
E y > 100 MeV (upper panel) and E y > 100 GeV (lower panel). 



dlogfi 



: R 3 p(R) 2 



R 3 R< R s , 
R- 3 R » R,. 



(43) 



Here we assumed a central plateau of the density profile which 
steepens beyond the scale radius R s and approaches the asymp- 
totic slope of R 3 of the dark matter profile tha t shapes the gas 
distribution at large radii dKomatsu & S eliak 2001, and references 
therein). This radial behavior makes the simulated y-ray luminosity 
only weakly dependent on uncertainties from the incomplete phys- 
ical modelling of feedback processes in the cluster cores. 

6.1 Contribution of different y-ray emission processes 

Figure [T5] shows the scaling relations of the IC and pion decay 
emission for two different energy scales of interest to the Fermi 
y-ray space telescope and imaging air Cerenkov telescopes. We 
compare the total emission and the contribution from the individ- 
ual emission components of each cluster. The very similar slopes of 
the mass-to-luminosity scaling relation at both energies (Table[5} is 
a consequence of the small variance in the proton spectrum (Fig[6) 
among different galaxy clusters. The individual emission processes 
also show similar slopes, with small scatter for the pion decay emis- 
sion and sIC component. Contrary, the pIC emission has a larger 
scatter than the secondary emission components due to the differ- 
ent dynamical states: the presence of strong merger or accretion 
shocks is critical for the generation of primary CR electrons and 
the associated radiative emission. 



Figure 16. Studying the galaxy bias to cluster y-ray scaling relations for 
energies E y > 100 MeV and R < R m . In the upper panel, the blue crosses 
show the cluster emission with galaxies and the red crosses show the clus- 
ter emission without galaxies, the solid lines show the mass to luminosity 
scaling (TablefS}- The lower panel shows the contribution from individual 
components (primary IC - red squares, secondary IC - green circles, 7T° 
decay - blue triangles) for emission including galaxies (light color) and ex- 
cluding galaxies (dark color). 



The ratio of the pion decay to the pIC emission in the 100 MeV 
and 100 GeV regime is very similar. This is partly a coincidence 
and owed to our particular choice of the two energy bands: the ef- 
fective spectral index of the pion decay between these two ener- 
gies is flattened due the pion bump. It happens to be similar to the 
power-law index of the pIC component which itself is unaffected 
by the energy cutoff of the electron spectrum at these energies. If 
we chose a smaller (larger) value than 100 MeV for the lower en- 
ergy band, we would obtain a lower (higher) pion-to-pIC ratio due 
to the steeper intrinsic spectrum of the CR protons at lower energies 
(cf.Hg.ID. 



6.2 y-ray emission from individual galaxies 

We investigate the bias of galaxies to the scaling relation of the y- 
ray luminosity above 100 MeV in Fig. [T6] The top panel shows 
the total y-ray emission, where the presence of galaxies biases 
smaller mass clusters slightly more compared to their larger ana- 
logues, with an average bias of about a factor two across our cluster 
sample. Masking galaxies reduces the overall scatter in the scaling 
relations, particularly at low masses. In the lower panel, we show 
the contribution from individual emission components. The largest 
bias originates from the pion decay (blue triangles) and the sIC 
emission (green circles), while the pIC component (red diamonds) 
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Figure 17. Influence of the emission from the accretion region around clus- 
ters on the mass to luminosity cluster scaling relations for E y > 100 MeV. 
In the upper panel, we show the luminosity integrated out to 3R v ; r (red 
crosses) as well as the luminosity integrated out to R yn (blue crosses). The 
solid lines indicate the mass to luminosity scaling relation found in TablefS] 
The lower panel shows the contribution from individual components (pri- 
mary - red squares, secondary - green circles, n° decay - blue triangles) for 
luminosity integrated out to R vu (light color) and luminosity integrated out 
to 3i? v i, (dark color). 

is barely affected by this masking procedure since it does not scale 
with density. 



6.3 y-ray emission from the cluster periphery and WHIM 

In this section we study the dependence of the accretion region 
around clusters on y-ray scaling relations. Specifically, we com- 
pare the total y-ray emission within R V1I to the emission within 
3 Ryi r which hosts the WHIM and individual satellite galaxies (and 
groups) that have not yet accreted onto the cluster. 

From Fig. [IT] it is clear that the total luminosity is dominated 
from the region inside R v [ r . For small clusters, the flux-correction 
from the WHIM is of the order of 30 percent which the correc- 
tion is smaller for massive systems - below 10 percent. This stems 
from the fact that the emission of low mass systems with smaller 
potential wells is easier to perturb - through accreting clumps of 
matter or nearby satellite systems. The pIC component contributes 
a factor 2-10 more in the WHIM than within R YII . The reason 
being that especially for merging systems, the pIC profile is rather 
flat (see Fig. Hence it contributes substantially to WHIM lu- 
minosity, whereas the density dependent secondary components in 
the WHIM are negligible due to the low gas densities. Only satellite 
systems within the WHIM contribute at a low level to the secondary 
components. This effect is especially pronounced when galaxies are 



excluded and implies that the pIC component becomes comparable 
to the pion decay emission for a few low mass clusters. Note that in 
the TeV regime, the pIC is considerably suppressed due to the lim- 
ited maximum energy of the primary CR electrons which reduces 
that effect in the WHIM. 



7 PREDICTION OF THE y-RAY EMISSION FROM 
NEARBY GALAXY CLUSTERS 

We use the mass-to-luminosity scaling relations as derived in Sec- 
tion [6] in combinations with the virial m asses of galaxy clus- 
ters o f the extended HIFLUGCS catalogue dReiprich & Bohringetl 
120021) - the "Highest X-ray FLUx Galaxy Cluster Sample" from 
the ROSAT all-sky survey - to predict their y-ray emission^ In 
Fig. ff8] we show the y-ray flux for radii R < R VK and energies 
above 100 MeV and 100 GeV, as a function of the identifier (ID) 
in the extended HIFLUGCS catalogue. We specifically name those 
clusters that have a flux (in our optimistic model) which is larger 
than the sensitivity of the Fermi all-sky survey after two years of 
data taking. We find that the brightest clusters in y-rays are Virgo, 
Ophiuchus, Coma, Perseus and Fornax. This result agrees well with 
previou s y-ray studies u s ing galaxy clusters from the HIFLUGCS 
sample jPfrommeJ[2 008; Jelte ma et alj2009l) . 

Figure[l8]should serve as a starting point to identify promising 
sources for the y-ray experiment in question. In a realistic setting, 
we would have to include the instrumental response and the point- 
spread function to obtain the predicted detection significance for 
the model in question. We note that this procedure of using scaling 
relations does not take into account deviations of individual sys- 
tems from the mean y-ray flux at a given cluster mass. One would 
rather have to model each system separately along the lines pre- 
sented in Section|5] 

To address the effect of source extension on the detection sig- 
nificance of Fermi, we compute the y-ray flux of each cluster within 
the Fermi equivalent angular resolution at 100 MeV in Fig. [19] To 
this end, we interpolate the scaling relations in Table [5] to the ra- 
dius corresponding to the angular resolution of 3.5 deg. We limit 
the size of each source to 3 J? v j r since there is negligible additional 
flux beyond this radius. For most clusters the flux is very similar to 
what we found in Fig.ff8]because of the similar mass to luminosity 
scaling relations for R vil and 3 R^. 



8 DISCUSSION AND COMPARISON TO PREVIOUS 
WORK 

8.1 Comparison to previous work on y-ray emission from 
clusters 

In support of the new instrumental capabilities in y-ray astronomy, 
several pioneering papers have appeared that simulate the high- 
energy y-ray emission from clusters. Here we make a comparison 
to some of those papers. 

Comparison to Pfrommer et al. - In t he ser ies of papers 
IPfrommer etal] d2007l l2008h and IPfrommen 1 j2008h simulate the 
same cluster sample as we do. In fact, our work represents an ex- 
tension of these earlier works. Overall, our results are in continuity 
with their results. The differences emerge from the details of the 



14 We have also added the Virgo cluster to the sample that we nevertheless 
refer to as extended HIFLUGCS catalogue in the following. 
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Table 5. Cluster y-ray scaling relations' 1 '. 



Notes: 





y-rays (E y 


> 100 MeV): 


y-rays (E y 


> 1 GeV): 


y-rays (E y > 


100 GeV): 


model* 2 ' 


r O) 


fh 


r <4) 


Py 


AS) 


Py 


Including galaxies, 1 R vlr 
Including galaxies, 2 i? v j, 
Including galaxies, 3 R v i, 
Excluding galaxies, 1 i? v j r 
Excluding galaxies, 2 R VII 
Excluding galaxies, 3 i? v j r 


5.24 ±0.99 
5.73 ± 0.97 
6.53 ± 1.20 
3.13 ±0.37 
3.32 ± 0.54 
3.44 ± 0.43 


1.34 ±0.09 
1.24 ±0.09 
1.14 ±0.08 
1.46 ±0.06 
1.41 ±0.08 
1.39 ±0.06 


7.83 ± 1.70 
8.48 ± 1.64 
9.63 ± 1.87 
4.02 ± 0.42 
4.22 ±0.53 
4.37 ± 0.57 


1.35 ±0.10 
1.24 ±0.10 
1.15 ±0.09 
1.43 ±0.07 
1.39 ±0.06 
1.37 ±0.07 


2.78 ±0.71 
3.00 ±0.71 
3.47 ± 0.76 
1.16 ±0.11 
1.24 ±0.19 
1.30 ±0.20 


1.33 ±0.13 
1.23 ±0.12 
1.14 ± 0.11 

1.34 ±0.08 
1.31 ±0.08 
1.29 ±0.08 



(1) The cluster y-ray scaling relations are denned by A = A M 1 ^, where M 15 = M vir /(10 15 M o ). 

(2) The y-ray luminosity from respective cluster is obtained by integrating over the region within an integer times R v j,. For scaling relations that include 
galaxies, we apply a central core cut-out and exclude a spherical region with radius r < 0.025 R v k that is centered on the brightest central y-ray point-source. 

(3) The normalization of the y-ray scaling relations (E y > 100 MeV) is given in units of 10 45 ph s , 

(4) The normalization of the y-ray scaling relations (E y > 1 GeV) is given in units of 10 44 ph s . 

(5) The normalization of the y-ray scaling relations (E y > 100 GeV) is given in units of 10 42 ph s . 
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Figure 18. Predicted y-ray flux in clusters and groups in the extended HIFLUGCS catalog to which we also add the Virgo cluster. For each cluster and 
group, we account for the flux from the region within one virial radius. The left axis shows the flux above 100 MeV while the right axis accounts for 
the flux above 100 GeV. The black line refers to our optimistic model where we include the flux contribution from galaxies and the red line shows the 
flux without galaxies (cf. Table|3J. We name the clusters and groups with T 7I o_ y (E y > 100 MeV) > 2 X 10 ph cirT 2 s in our optimistic model which 
roughly corresponds to the sensitivity of the Fermi all-sky survey after two years of data taking. 



CR physics, where they adopted a simplified description using a 
single CR population with a spectral index of 2.3. Hence they con- 
centrated on the non-thermal radio emission as well as the y-ray 
emission at energies E y > 100 MeV that depend only weakly on 
the particular value or even a running of the CR spectral index (as 
long as it is close to the true one). For the primary electron popu- 
lations, a maximum electron injection efficiency of £ e = 0.05 was 
assumed. In addition, the bias from anomalous galaxies was not 
addressed. 

Comparing profiles of the total surface brightness above 
100 MeV in our optimistic model w ith galaxies, shown in grey 
in Fig. [4] to the brightness profiles in IPfrommer et alj d2008l) . we 
find only very small differences. These differences are caused by 
a combination of different binning and Poisson noise in the galax- 
ies' spatial distribution which have a different realization due to our 



slightly modified CR description that changes the hydrodynamics. 
The largest difference is seen for the pIC component in the periph- 
ery of merging clusters and a consequence of the different values 
for the injection efficiency f e adopted. In addition, we compare the 
mass to luminosity scaling relation above 100 MeV in our opti- 
mistic model (see Table [5} a nd find tha t they agree to the percent 
level with what was found in Pfrommet] d2008l) . 

Comparison to Miniati et al. - There have been a series of pi- 
oneering papers simulating the non-thermal emission from clusters 
by numerically modelling discretized CR proton and electro n spec- 
tra on top of Eulerian grid- b ased cosmological si mulations dMiniatil 
l200lMMiniati et alj200lj|bllMiniatil2002ll2003l) . In contrast to our 
approach, these models neglected the hydrodynamic pressure of the 
CR component, were quite limited in their adaptive resolution ca- 
pability, and they neglected dissipative gas physics including radia- 
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Figure 19. Predicted y-ray flux above 100 MeV in clusters and groups in the extended HIFLUGCS catalog to which we also add the Virgo cluster. The 
flux comes from the region within the Fermi angular resolution at 100 MeV, i.e. a circular region of radius 3.5 degree that contains 68 per cent of the PSF, 
but with the limit at 3 fi v j r for each cluster and group. The black line refers to our optimistic model where we include the flux contribution from galaxies 
and the red line shows the flux without galaxies (cf. Table|5J. We name the clusters and groups with T^_ y (E y > 100 MeV) > 2 X lO^ph cm 2 s 1 in 
our optimistic model which roughly corresponds to the sensitivity of the Fermi all-sky survey after two years of data taking. 



tive cooling, star formation, and supernova feedback. Comparing 
the y-ray emission characteristics of the IC emission from primary 
CR electrons and hadronically generated secondary CR electrons as 
well as the pion decay y-rays, we confirm the qualitative picture of 
the emission characteristics of the different y-ray components put 
forward by these authors. However, we find important differences 
on smaller scales especially in cluster cores, the emission strength 
of the individual components and their spectra. 

We confirm that the high-energy y-ray emission (E y > 
100 MeV) from cluster cores is dominated by pion decays while 
at lo wer energies, the IC emission of secondary CR electrons takes 
over jMiniatil2003h . We reproduce their finding that the y-ray emis- 
sion in the virial regions of clusters and beyond in super-cluster 
regions is very inhomogeneous and stems in part from the IC emis- 
sion of primary shock accelerated electrons. Contrarily to these au- 
thors, we find that the surface brightness of this emission compo- 
nent remains sub-dominant in projection compared to the hadroni- 
cally induced emission components in the cluster core and that the 
pion decay completely dominates the high-energy y-ray emission 
of clusters above a few MeV (cf. Fig. QJ. In addition we predict 
a pIC spectrum that is somewhat steeper with a photon index of 
r =! 1.15 which resembles a steeper primary electron spectrum 
in our simulations compared to theirs. This points to on average 
weaker shocks that are responsible for the acceleration of primary 
CR electrons that dominate the pIC emission. This discrepancy of 
the pIC spectral index causes the discrepancy of the pIC flux at high 
y-ray energies. 

In the WHIM we find that the pIC emission dominates the total 
y-ray spectrum below about 100 MeV, and a comparable flux level 
of n° -decay and pIC between 100 MeV and 1 TeV, where the n°- 
decay takes over. This is in stark contrast to the finding of Miniati 
(2003), where the pIC is dominating the pion decay emission by a 
factor of about 10 or more over the entire y-ray energy band. We 
note that our y-ray fluxes from clusters are typically a factor of two 



smaller than the estimates given in iMiniatiet allfeOOld) which has 
important implications for the detectability of clusters by Fermi. 

There are several factors contributing to the mentioned dis- 
crepancies. (1) Our simulations are Lagrangian in nature and hence 
adaptively resolve denser structures with a peak resolut ion of 
5/T 1 kpc. In contrast, the cosmological simulations of Miniati 
(2003) have a fixed spatial resolution of ~ 100 ft -1 kpc which is 
too coarse to resolve the observationally accessible, dense cen- 
tral regions of clusters in this grid-based approach and underesti- 
mates CR cooling processes such as Coulomb and hadronic losses. 
It also c annot resolve t he adi abatic compression of CRs into the 
core. (2) iMiniati et alj d2000l) identified shocks with Mach num- 
bers in the range 4 <, M < 5 as the most i mportant in thermalizin g 
the plasma. In contrast. | Rvu et iD d2003h . |Pfrommer et alj J200r3) , 
and lSkillman et al] ( bOOa) found that the Mach number distribution 
peaks in the range 1 < M < 3. This finding seems to be robust 
as different computational methods have been used which range 
from fixed and adaptive Eulerian grid codes to Lagrangian Tree- 
SPH codes. Since diffusive shock acceleration of CRs depends sen- 
sitively on the Mach numb er, this implies a mor e efficient CR in- 
jection in the simulations bv lMiniati et alj d2001al) . It also results in 
a flatter CR electron and CR ion spectru m comp a red to ours shown 
in Fig. [T] Hence, the pIC emission of Miniati (2003) has a flat- 
ter photon ind ex and a boosted flux. (3) For the CR ion spectrum, 
lMiniatilj2003h uses only four momentum bins which is not enough 
to resolve the pion bump accurately. The large pion decay plateau 
which he found indicates a constant CR ion spectral index in this 
energy range. This is in contradiction to the concavely shaped CR 
spectrum that our cluster simulations show, where the shape is a 
consequence of the Mach number statistics and the adiabatic trans- 
port. The difference in the CR spectral shape is especially important 
for CR energies above 1 GeV, since those CRs give rise to the n°- 
decay emission at energies above the pion bump. 

Comparison to Kushnir et al. - iKushnir & Waxmaril d2009l) 
use a simple analytic model to follow the evolution of ICM CRs, 
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accelerated in strong accretion shocks. Interestingly, their approach 
predicts similar characteristics for the pion decay emission, in par- 
ticular its flux agrees with our prediction within a factor two. In 
contrast, their model predicts a high-energy y-ray flux of the pIC 
component that is approximately a factor of 600 larger than ours 
due to a different spectral description where they adopted a spectral 
index of 2. This, however, is in conflict with our simulated average 
spectral index of 2.3 that is a conseque nce of cluster assembly his - 
tory (compare our Fig.|4]with Fig. 3 in Kushnir & Waxman 2009). 
This finding, together with different adopted values for the elec- 
tron injection efficiency and CR-to-thermal ratio, led them to the 
contradicting conclusion that the expected overall y-ray emission 
would be much more extended. We believe that their model would 
over-predict the amount of observed radio relic emission in clus- 
ters, in particular when considering magnetic field amplification at 
accretion shocks on a level that is o nly a fraction of wha t is ob- 
served in supernova remnant shocks dUchivama et alj|2007l) . It can 
be easily seen that the different power-law indices are indeed the 
reason for the flux discrepancies by comparing the energy flux of 
electrons at 170 GeV which are responsible for IC photons at 100 
MeV (assuming the up-scattering of CMB photons). Adopting our 
effective spectral injection index of primary electrons of or c _j n j = 2.3 
and assuming a post-shock temperature at a typical accretion shock 
of kT =a 0. 1 keV, we find a flux ratio between their model and ours 
of (0. 1 keV/170 GeV) 03 ^ 600. 

We would like to compare their model predictions for the pion 
decay emission in more detail. To this end, we use a proton injec- 
tion efficiency in their model of rj p a 0.2 for the comparison. In 
our simulations we use a maximum proton injection efficiency of 
'Jmax.p = 0.5. The proton injection efficiency in our simulations is 
dynamical, and depends on the strength of the shocks. Since the 
cluster evolves with time and the majority of CRs are injected at 
higher redshift during the formation of clusters, the final CR pres- 
sure depends on the interesting interplay of the actual value of the 
shock injecti on efficiency and t he successive CR transport. Using 
results from fenfilin et alj2007h . we estimate an effective injection 
efficiency of approximately rj p = 0.2. We note, however, that their 
baseline model assumes 7/ p = 0.02, which will suppress the flux by 
a factor 10 in their model. We now contrast the y-ray flux predic- 
tions of the Coma and Pers eus cluster of our semi-anal ytic model 
to the one worked out in iKushnir & Waxmanl {2009). First we 
study the y-ray flux from Coma within R s = R yil above 100 MeV 
where we find with their model a flux 7v>_ yk&w (> 100 MeV) = 
(0.70 - 1.7) x 10~ 9 phcnr 2 s" 1 . This flux is only a factor two lower 
than what we predict for the pion decay emission from Coma. Turn- 
ing to pIC, which play a much more important role for the total 
y-ray flux than the 7r°-decay emission in their model, results in the 
flux r p ic,k&w(> 100 MeV) a 1.5 x 10~ 8 phcnr 2 s _1 . This is only a 
factor two below the EGRET upper limit and can be readily tested 
with the one-year data from Fermi. 

Studying fluxes in the TeV y-ray regime is also of great im- 
portance since it compares the spectral representation of t he mod- 
els. Using the analytic model o f lKushnir & Waxmart {20091) , results 
in the flux r^o_ ycomil (> 1 TeV) = 9.3 x 10" 14 phcnr 2 s" 1 within 
9 = 0.2 deg that is about a factor 10 below the upper limit for 
Coma set by HESS and only a factor two larger than the flux we 
predict for Coma. However, their result is most probably flawed by 
their too simplistic assumption for the CR spectral index of a = 2. 
If we use our universal concave shaped spectrum instead of their 
flat CR spectrum, we can show that their flux would decrease by a 
factor ~ (1 TeV/100MeV) <r ™Mev-" ~ io 4 * - 24 ~ 9. 

Finally we study the surface brightness profiles from Coma 



and Perseus predicted by our semi-analytic model and compare 
it to theirs. The dashed and dotted-dashed lines in Fig. [13] show 
their predictions for the Coma cluster above the energy 100 MeV 
and 100 GeV, respectively. Using their formalism we find a surface 
brightness above 100 MeV that is about a factor two smaller than 
what we predict. However, above 100 GeV our predictions are in 
better agreement. 

8.2 Limitations and future work 

The ideal CR formalism would trace the spectral energy evolu- 
tion, as well as the spatial evolution of CRs, and at the same time 
time keep track of the dynamical non-liner coupling with magneto- 
hydrodynamics. In order to make cosmological simulations less ex- 
pensive in computational power, we are forced to make compro- 
mises. The simplifying assumptions chosen, enable us to run cos- 
mological simulations of the formation of galaxy clusters with the 
necessary resolution to resolve their cores. At the same time, these 
assumptions enable us to follow the CR physics self-consistently 
on top of the radiative gas physics. Here we outline our most se- 
vere limitations for computing the y-ray emission from clusters. 

(i) In our simulations we neglect the effect of microscopic CR 
diffusion and CR streaming. The collisionless plasma forces CRs to 
stay predominantly on a given field line and to diffuse along it. The 
random walk of field lines cause initially closely confined CRs to 
be transported to larger scales which can be described as a diffusion 
process. In our model we assume the magnetic field to be tangled 
on scales smaller than those we are interested in, A ~ 10 kpc in the 
center and even larger scales outside. Hence, CRs are magnetically 
coupled to the thermal gas and advected alongside it. The diffusiv- 
ity can be rewritten into a macroscopic advection term that we fully 
resolve in our Lagrangian SPH simulations by construction and a 
microscopic diffusivity. The advection term dominates over micro- 
scopic term, as the following estimate for the diffusivities shows: 
K^ y = 100 kpc x 1000 km/s 10 3U cm 2 /s » k as 10 29 cm 2 /s. 
Further work is needed to study microscopic anisotropic diffu- 
sion, in combination with self-consistent modelling of the magnetic 
fields. 

(ii) We also did not account for the injection of CRs by AGN or 
supernova remnants where the additional CRs would diffuses out of 
AGN-inflated bubbles or drive starburst winds that enrich the IGM. 
In addition we do not account for the feedback processes by AGN 
despite their importance for understanding the nature of the very 
X-ray luminous cool cores found in many clusters of ga laxies. For 
further details we refer the reader to lSiiacki et al.[ (2008). 

(iii) We postpone the study of the potential contribution of a 
population of re-accelerated electrons to the IC y-ray emission 
throughout this work: strong merger shocks and shear motions at 
the cluster periphery might inject hydrodynamic turbulence that 
cascades to smaller scales, feeds the MHD turbulence and even- 
tually might be able to re-accelerate an aged CR electron popu- 
lation. Due to non-locality and intermittency of turbulence, this 
could partly smooth the very inhomogeneous primary emission 
component predominantly in the virial regions of clusters where 
simulations indicate a higher energy density in random motions. 
However, to study these effects, high-resolution AMR simulations 
are required that refine not only on the mass but also on some 
trac er for turbulence such as the dimensionless vorticity parame- 
ter {iapichino et alj200 8: Iapichino & Niemever 2008; Maie r et all 
120091) . 

(iv) Our model for the diffusive shock acceleration assumes a 
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featureless power-law for both, the proton and the electron acceler- 
ation, that is injected from the thermal distribution. The complete 
theoretical understanding of this mechanism is currently an active 
research topic that includes non-linear effects and magnetic field 
amplification dVladimirov et alj l2006h . Phenomenologically, we 
believe that there are strong indications for the diffusive shock ac- 
celeration mechanism to be at work which come from observations 
of supernova remnants over a wide ran ge of wavelengths from the 
radio , X-rays into the TeV y - ravs (e.g..lEllisonll20 00; Hughes et al. 



200d; lEllison & et al.1 12003; IWarren et alj 120051 : lAharonian et al. | 
2004 , 120061) as well as the bow shock of the Earth (Ellison et al. 



199C ; Ishimada et a il l 19991) . Theoretical work suggests that the 



spectrum of CRs which is injected at strong shocks shows an in- 
trinsic concave curvature: the feedback of the freshly accelerated 
and dynamically important CR pressure to the shock structure re- 
sults in a weaker sub-shock that is proceeded by a smooth CR 
precursor extending into the upstream. Hence low-energy protons 
are only shock-compressed at the weaker sub-shock and experi- 
ence a smaller density jump which results in a steeper low-energy 
spectrum (compared to the canonical value a = 2 from linear the- 
ory). In contrast, the Larmor radii of high-energy protons also sam- 
ple the CR precursor and experience a much larger density con- 
trast that results in a flatter high-energy spectrum with a < 2 
jEllison &Eichlerl 1 19841 ; lAmato et alj|2008l : ICaprioli et ai]l2009l) . 
The low-energy part of the CR spectrum in clusters (as found in 
this work) should be unaffected since a softer population of CRs 
dominate there with a ~ 2.5 and non-linear effects are presum- 
ably negligible in this regime. However the high-energy part of the 
CR spectrum in clusters could become harder compared to what we 
found due to these non-linear effects. Future work will be dedicated 
to improve our model and to incorporate more elaborate plasma 
physical models and to study the uncertainty of our results with re- 
spect to the saturated value of our CR a cceleration efficiency (e.g., 
iKang & Jonesll2007l; lEdmon et alj|2007l) . 

(v) An artificial surface tension effect limits the ability of SPH 
(in its standard conservative form) to follow the growth of bound- 
ary instabilities such as the Kelvin-Helmholtz instability at the in- 
terface of a dense and unde r-dense phase accu rately, i.e. on the 
predicted linear growth time dAgertz et alj|2007l) . In the context of 
galaxy cluster simulations, this only occurs at the interface of the 
ISM of individual galaxies and the ICM. This causes an unphysi- 
cally long survival time of dense gaseous point sources after they 
got ram pressure/tidally stripped from their galactic halo - for sim- 
plicity, we call them "galaxies" and describe the physics in detail in 
Section |231 In our paper, we decided to show our result for an opti- 
mistic model that includes all galaxies and one conservative model 
that cuts all galaxies. This is meant to bracket the realistic case. 
Also, the main result is based on our conservative model where we 
cut out the galaxies. In this way we circumvent these issues. 

8.3 Impact of these limitations on our results 

Generally, we acknowledge that the CR spatial distribution is more 
uncertain than the spectrum due to the details outlined in Sec- 
tion [jO] Below, we detail our considerations why we believe that 
the spectrum that we found is robust even when considering uncer- 
tainties such as additional CRs injected from AGN, re-acceleration 
of CRs at MHD turbulence, CR diffusion, and non-linear shock ac- 
celeration. We outline here the reasons in detail: 

Additional CRs injected from AGN. It is very uncertain 
whether AGN j ets are powered hadronically or though Poynt- 
ing flux (e.g., ICelotti et alj 1 19981: Iffirotani & Okamotol 1 19981 : 



Sikora & Madejski 2000). Irrespectively, the energetics of AGN are 
insufficient to account for a majority of CRs in clusters - in particu- 
lar for the most massive systems (Thompson & Pfrommer in prep.). 

CR re-acceleration through MHD turbulence. The involved 
physics is currently very uncertain such as the level and na- 
ture of turbulence in the ICM, how CRs exactly interact with 
plasma waves, how efficient this accelerates CRs, and whether a 
power-law extrapolation between the gyro radius of a CR, R ~ 
10~ 5 pc (B/yuG) -1 (.E/lOGeV), and the scales accessible to current 
simulations with peak resolution of a few kpc is justified. Hence 
it appears that is is impossible to constrain the impact of turbu- 
lent re-acceleration on the CR spectrum in clusters self-consistently 
from first principles. However, in our Milky Way, we are able to 
understand the observed CR spectrum on Earth with /cr ~ p~ 2 7 
fairly well in terms of injection and transport. CRs with an injected 
spectrum of p~ a 3-2 - 4) experience momentum dependent diffusion 
so that the more energetic particles can leave the system in a so- 
called 'leak y-box model'; an e ffect that accounts for the observed 
steepening lschlickeisej|2002h . This leaves little room for spectral 
modifications through turbulent re-acceleration. Hypothesizing that 
the fundamental interactions of plasma waves with particles should 
not be very different in the ISM and ICM, we believe that we are 
safe to neglect this process to first order. Note, however, that this 
argument neglects possible important CR transport processes that 
might become important in the cluster environment due to the much 
longer CR life time compared to the ISM in our Galaxy. 

Spatial diffusion of CRs. If spatial CR diffusion is momen- 
tum dependent it will introduce a radial dependence in the shape of 
the CR spectrum since high energy CRs can diffuse out of the cen- 
tral regions faster. This, in turn, will affect the observed morphol- 
ogy and spectrum of all relevant non-thermal emission components, 
including gamma-rays from pion decay, IC from secondaries, and 
synchrotron emission. 

Non-linear shock acceleration. Since the CR spectrum at 
GeV-to-TeV energies is sufficiently steep, intermediate Mach num- 
ber shocks are responsible for the acceleration of these CRs - with 
efficiencies there are modest (not in the saturated regime) so that the 
non-linear back-reaction is expected to be small or even negligible. 

Over-cooling problem. Due to the over-cooling problem, the 
modeling of cluster cores with radiative cooling but without any 
significant feedback process produces an unphysically high stel- 
lar mass fraction. This is not expected to impact the CR spectrum 
significantly. However, the effect on the CR morphology might be 
substantial and might depend on the details of the required feed- 
back process. This can potentially impact the non-thermal cluster 
observables. Future work is needed to solve this problem. 



9 CONCLUSIONS 

In this paper we have simulated 14 galaxy clusters spanning two or- 
ders of magnitude in mass and a broad range of dynamical stages. 
The simulations follow self-consistent CR physics on top of the 
dissipative gas physics including radiative cooling and star forma- 
tion. We have simulated high-energy y-ray emission maps, profiles 
and spectra of various emission components. These include the in- 
verse Compton emission from primary, shock-accelerated electrons 
(pIC) and secondary electrons that result from hadronic interactions 
of CR protons with ambient gas protons (sIC), as well as y-rays 
from neutral pion decay that are also generated in these hadronic 
reactions. 

We would like to emphasize that we focus on the intrinsic 
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spectrum emitted at the cluster position without taking into ac- 
count photon propagation effects to highlight the various physical 
process that shape the emission spectra. Depending on the cluster 
redshift, these spectra attain a high-energy cutoff due to e + e~-pair 
production on IR and optical photon s which can be easily de rived 
from the photon-photon opacity (e.g.. lFranceschini et al.l2008l) . We 
also caution the reader that we assume an optimistic value for the 
maximum shock injecti on efficiency (base d on data from super- 
nova remnant studies bv lHelder et alj|2009h ; smaller values would 
reduce the resulting y-ray emission accordingly. To date it is not 
clear whether these high efficiencies apply in an average sense to 
strong collisionless shocks or whether they are realized for struc- 
ture formation shocks at higher redshifts. Hence the goal of this 
work is to establish a thorough framework and to predict the level 
of y-ray emission that we expect for this efficiency. We note that 
one cannot lower the acceleration efficiency infinitely if one wants 
to explain radio (mini-)halos in the hadronic model of CR interac- 
tions. For clusters that host such a large, unpolarized, and centrally 
peaked radio halo emission that resembles the thermal X-ray sur- 
face brightness, one can derive a minimum y-ray flux. The idea is 
based on the fact that a steady state distribution of CR electrons 
loses all its energy to synchrotron radiation for strong magnetic 
fields (B » Bcmb 3.2;uG x (1 + z) 2 ) so that the ratio of y-ray to 
synchrotron flux becomes independent of the spatial distribution of 
CRs and thermal gas. Lowering the magnetic field would require 
an increase in the energy density of CR electrons to reproduce the 
observed synchrotron luminosity and thus increase the associated 
y-ray flux (for applicatio ns to Coma and Perseus, see IPfrommeJ 
20Q8s: lAleksic et alfeoid respectively). 

According to our simulations, clusters have very similar mor- 
phology in the 100 MeV - 100 GeV Fermi band, and in the 
100 GeV - 10 TeV Cerenkov band. This is due to the power- 
law spectra of the dominating pion decay emission (which show 
a slowly running spectral index) and ultimately inherited by the 
parent CR proton distribution. The emission from the central parts 
of clusters are dominated by y-rays from pion decay, while the pe- 
riphery of the ICM and the WHIM have a considerable contribution 
from pIC, which is especially pronounced in merging clusters. The 
energy dependent photon index for 100 MeV to 1 GeV has a me- 
dian value of T = 0.9 due to pion decay induced emission in the 
central parts of the clusters, while that in the periphery shows a 
slightly higher value of T = 1.1 which is due to the substantial con- 
tribution from pIC. In the energy range from 100 GeV to 1 TeV, 
the photon index steepens to T = 1.25 in the central regions. This 
spectral steepening in the cluster center is due to the convex cur- 
vature of the pion bump around 100 MeV causing a steepening in 
the asymptotic y-ray spectrum at higher energies. The small con- 
cave curvature at higher energies is not able to compensate for this 
effect. At energies E y > 1 TeV, the photon index in the cluster out- 
skirts attain a much higher value due to a super-exponential cutoff 
of the primary IC spectrum. This emission component contributes 
substantially to the total y-ray emission there. At these energies, 
the electron cooling time is smaller than the time scale for dif- 
fusive shock acceleration which causes this cutoff in the electron 
spectrum. We used a semi-analytic formula for the injected elec- 
trons from which we derive the cooled electron distribution with 
the characteristic super-exponential cutoff. The shape of this spec- 
trum is passed on to the pIC spectrum and we capture this shape 
with a fit that is valid both in the low-energy Thomson regime as 
well as in the high-energy Klein-Nishina regime. 

The simulated CR proton spectra show an approximate power- 
law in momentum with a few additional features; a cutoff at p = 



P p /m p c 2 =! 0.1, a concave shape between p ^ 1 — 10 6 , and a steep- 
ening by p 1/3 between p 10 6 - 10 s . The overall shape of the 
spectrum shows only little variance between the clusters, indicat- 
ing a universal CR spectrum of galaxy clusters. The radial depen- 
dence of the spectrum within the virial radius is negligible to first 
order. This allowed us to construct a semi-analytic model of the 
median CR proton spectrum across our cluster sample. Using the 
semi-analytic CR spectrum we derive a semi-analytic formula for 
the y-ray flux from the pion decay induced emission that dominate 
the total y-ray spectrum above 100 MeV. We apply this formalism 
to the Perseus and Coma clusters, using their density profiles as 
inferred from X-ray measurements and predict that the flux from 
Perseus is clo se to the recent uppe r limits obtained by the MAGIC 
collaboration jAleksic et alj201Ch . 



The mass-to-luminosity scaling for the 100 MeV, 1 GeV, and 
100 GeV regimes show very similar slopes for both the total y- 
ray luminosity and all the components, which is due to the small 
variance in the CR spectrum. Masking galaxies decreases the to- 
tal y-ray emission by a factor of 2-3. The cut has a larger effect 
on smaller mass clusters since the emission of low mass systems 
with smaller potential wells are easier to perturb - through accret- 
ing clumps of matter or nearby satellite systems. We also found 
that the presence of galaxies considerably increases the scatter in 
the y-ray scaling relation. The region outside R V1I only contributes 
marginally (of order ten per cent) to the total y-ray emission for 
massive clusters while it contributes significantly to the total y-ray 
luminosity of low mass clusters with a factor < 1.5. This is again 
mostly due to the pion decay emission from satellite systems that 
have not yet accreted on the cluster. The flux of the pIC compo- 
nent is increased by a factor of 2 - 10 when the WHIM is included. 
This can be explained by the rather flat spatial profiles of the pIC 
emission. 



Combining our y-ray scaling relations with the virial masses 
of galaxy clusters of the extended HIFLUGCS catalogue, we pre- 
dict a detection of a few galaxy clusters above 100 MeV with Fermi 
after two years, where Virgo, Ophiuchus, Coma, Perseus and For- 
nax are expected to be the brightest clusters in y-rays (barring un- 
certainties in the injection efficiency). Since Fermi already discov- 
ered the central AGN in Virgo/M87 and Perseus/NGC1275 the de- 
tection of the somewhat more extended and dimmer pion decay 
component will be very challenging in these clusters and requires 
careful variability studies to subtract the AGN component. For en- 
ergies above 100 GeV, the flux of these clusters as determined by 
our scaling relation is more than 5 x 10 3 times lower. This provides 
a challenge for current Cerenkov telescopes as it is almost an order 
magnitude lower than the 50 h sensitivities. However, future up- 
grades of IACTs or the CTA telescope might considerably change 
the expectations. We note however that these estimates are too con- 
servative for cool core clusters, which are known to show enhanced 
X-ray fluxes by a factor of up to ten relative to clusters on the X-ray 
luminosity scaling relation. Since we expect the X-ray luminosity 
to tightly correlate with the y-ray luminosity, this sub-class of clus- 
ters should provide very rewarding targets due to the ample target 
matter for inelastic collisions of relativistic protons leading to y- 
rays. Applying our semi-analytic model for the y-ray emission, we 
identify Perseus among the best suited clusters to target for the cur- 
rent IACT experiments. 
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APPENDIX A: SUPPORTING MATERIAL FOR OUR 
SEMI- ANALYTIC MODEL: SPECTRUM AND SPATIAL 
DISTRIBUTION 

Here we present additional details on the spectral and spatial distri- 
bution of cosmic rays (CRs) that are important for the consistency 
of our semi-analytic model. 

Al Formalism for the semi-analytic modelling 

The CR spectrum for each SPH particle at the dimensionless proton 
momentum p = P p /m p c = 1 is given by f(p = l,R) = C(R), 
where we assume that the low-momentum cutoff q < 1. The CR 
normalization is denoted by C = Cp/m p , where C has units of 
inverse volume. At each radial bin, we use the the volume weighted 
C to calculate the normalized spectrum through 
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where the sum extends over SPH particles labeled by ('. In our semi- 
analytic formalism, we provide a fit to the mass weighted C, de- 
noted by Cm, and use the gas density profile (which in fact is the 
volume weighted density p v , that we simply denoted by p through- 
out the paper). The two methods to calculate f y are equivalent for 
all radial bins, since 
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We use the CR proton spectrum to calculate the integrated y- 
ray source density A by integrating the inverse Compton (IC) and 
pion decay y-ray source functions s y (E y ) giv en by equation (43) 
and by equation (19) of IPfrommer & Enfflinl J2004I) . respectively. 
The luminosity is calculated through the volume integral of A. In 
analogue to the CR spectrum, we show here that the pion decay 
induced luminosity on an SPH basis for any radial bin, 
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In the semi-analytic expression in equation JA4ft . we have used that 
the spectral part is separable from the spatial part in the first approx- 
imation, and that C is only a weak function of radius in the second 
approximation. Both these assumptions are validated by the uni- 
versal CR spectrum and approximate spatial universal profile that 
we found across our cluster sample where the details are shown in 
Section|4] 
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Figure Al. Convergence test for the CR spectrum: the main panel shows 
the volume weighted CR spectrum as a function of dimensionless CR mo- 
mentum p = Pp/m„c for the small CC cluster g676. The red line rep- 
resents our standard spectral resolution using a = (2.1,2.3,2.5,2.7,2.9), 
and the blue line shows the finer spectral resolution given by a = 
(2.1, 2.2, 2.3, 2.4, 2.5, 2.7). The bottom panel shows the difference between 
the two simulations normalized with the finer spectral resolution simulation. 



A2 Spectral convergence test 

Here we test the convergence of our updated CR model that we 
use in our simulations, where we account for both the spatial and 
spectral information of the CRs. In particular we allow for multiple 
CR populations, where each population is characterized by a fixed 
spectral index. For the test we run simulations with a finer spacing 
between the spectral indices, a = (2.1,2.2, 2.3,2.4,2.5,2.7), where 
we omitted the last bin a = 2.9 that was found to have negligible 
contribution to the spectrum. In Fig. lAll we show the result of the 
convergence test of the CR spectrum for the small CC cluster g676. 
The red line shows g676 with our standard spectral resolution and 
the blue line shows the g676 simulation with a finer spectral resolu- 
tion. The difference between the two simulations normalized with 
the spectrum from the better spectral resolution simulation is shown 
in the lower panel. In the GeV and TeV region the difference is less 
than 10 percent, showing that we are accurately able to capture the 
CR spectrum with our choice of a that has a wider spacing between 
the spectral indices. 



A3 Variance of the spatial CR distribution 

In this section we discuss the details of our spatial semi-analytic 
modelling and address the particle-by-particle variance of the spa- 
tial CR profile within a cluster. To this end, we plot the correlation 
space density of the dimensionless CR normalization C versus the 
overdensity of gas (5 gas = p/(Qbp CI ) - 1 for individual clusters. We 
over-plot the mean together with the 1-sigma standard deviation 
Cm as a function of overdensity. The scatter of 0.3 dex is roughly 
constant with density. This variance is most probably caused by 
the variance in shock strength and associated CR acceleration effi- 
ciency among different fluid elements in the past history of a clus- 
ter. If CRs are adiabatically transported into a cluster, we expect 
a weak radial dependence of the dimensionless normalization of 
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Figure A2. Correlation space density between the dimensionless normal- 
ization of the CR spectrum, C, and the overdensity <5 gas = p/(Slbp cr ) - 1, 
both in logarithmic representation. Shown is our simulated post-merging 
g72a cluster. The solid line shows the mean Cm as a function of the over- 
density and the error bars the 1-sigma standard deviation. 



the CR spectrum, C ~ p (a -»' 3 dEnfflin et al.ll2007h . However, we 
stress that the scale dependence of DM halos and which seems to 
be inherited by the gas and CR distribution could have an impor- 
tant effect in shaping C. Similarly, non-adiabatic CR transport ef- 
fects could have another important impact on C. The outer cluster 
regions (particularly of mergers) typically host weak to intermedi- 
ate strength shocks that increase C in that region. In contrast, in the 
very central cluster regions, the hadronic losses dominate and cause 
a suppression of C. 



APPENDIX B: PRIMARY, SECONDARY ELECTRONS, 
AND INVERSE COMPTON EMISSION 

In Section IBTI we introduce the steady state cosmic ray electron 
(CRe) population for secondary electrons. These electrons are cre- 
ated through CR proton-proton interactions via charged pions de- 
caying. The primary electrons in our simulations are accelerated 
through diffusive shock acce leration using the thermal leakage 
model originally proposed bv lEllison et alj ( Tl98ll) . In Section lB2l 
we start by deriving the primary steady state CRe population under 
the assumption that the diffusive time scale can be neglected. Then 
we continue the discussion for primary electrons where we addi- 
tionally account for an energy dependent diffusion that in combina- 
tion with the losses leads to a cutoff in the electron spectrum. In lB3l 
we use these CRe spectra to calculate the integrated secondary and 
pr imary IC source fun ctions. This description follows the approach 
of EnBlin et al. ( 2007), but we refer the reader to the appendix of 
Pfrommer et al. Il200j) for a detailed discussion about electron ac- 
celeration, cooling, and IC emission within the framework of the 
SPH formalism. 



Bl Secondary electrons 

Here we provide the equilibrium distribution of secondary CR elec- 
trons above a GeV. It is shaped by the high-energy part of the CR 



proton population that interacts with ambient gas and creates elec- 
trons through hadronic reactions. These electrons cool through IC 
and synchrotron radia tive losses. The CRe equilibrium spectrum 
jPfrommer et al J[2008h is derived through the balance of injection 
and losses, and is steeper by one compared to the CR proton spec- 
trum (a e j = a, r + 1). The volume weighted spectrum is given by 
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where Cm(R) is the mass weighted CR proton normalization (equa- 
tion I22t . A, the relative normalization for each CR population 
(equation|26t, and cr pp ; denote the effective inelastic cross-section 



for proton-proton interactions and is defined in Section pTl In addi- 
tion, o"t represents the Thompson cross-section, e b is the magnetic 
energy density, and £ ph denotes the photon energy density, taken to 
be that of CMB photons. 



B2 Spectrum of shock-accelerated electrons 

In clusters of galaxies, the dynamical and diffusive time scales of 
electrons are much longer compared to the shock injection and 
IC/synchrotron time scales. The radio synchrotron emitting elec- 
tron population cools on such a short time scale T sync < 10 s yrs 
(compared to the very long dynamical time scale T dyn ~ 1 Gyr) 
that we can describe this by instantaneous cooling at each timestep 
- in contrast to the CR protons. In combination with the fact that 
the CRes have a negligible pressure contribution, this enables us to 
account for the CRes in the post-processing. In this instantaneous 
cooling approximation, there is no steady-state electron population 
and we would have to convert the energy from the electrons to in- 
verse Compton and synchrotron radiation. Instead, we introduce a 
virtual electron population that lives in the SPH-broadened shock 
volume only, defined to be the volume of energy dissipation. Within 
this volume, which is co-moving with the shock, we can use the 
steady-state solution for the distribution function of CR electrons 
and we assume no CR electrons in the post-shock volume, where 
no energy dissipation occurs. Thus, the CR electron equilibrium 
spectrum can be derived from balancing the shock injection with 
the IC/synchrotron cooling: above a GeV and below 30 TeV, it is 
given by 



/ c (p e )dp c 
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Here, we introduced the unit-less electron momentum p e = 
P c l(m c c), where P c is the electron momentum. The spectral index 
of the equilibrium electron spectrum is denoted by a e = or™ + 1, 
where is the spectral index of the injected electron population in 
one-dimensional momentum space given by equation (0. The CR 
electron normalization scales linearly with the gas density C c oc p 
which we evolve dynamically in our simulations and depends in- 
directly on a in j and the dissipated energy rate per electron, E iiss , 
through the normalization of the injected CRes, C m ^ (for further de- 
tails see lPfrommer et ai1l2008l) . Tj„j represent the electron injection 
time scale, which depends on the time it takes for an electron to 
pass through the broadened shock. 

The shape of the steady-state electron power-law spectrum in 
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equation changes when the energy of the accelerated elec- 
trons reach a maximum electron energy that is determined by the 
competition between the diffusive acceleration, and radiative syn- 
chrotron and IC losses. In the cutoff region, the spectrum is pro- 
portional to the product of two terms - a power-law term (p c '~ 
which includes the phase space volume) and a super-exponential 
term, exp(-pl/ p 2 wdx ). The first term reflects a pile-up of electrons 
as their cooling time becomes comparable to the acceleration time. 
The exponential term, however, effectively cancels this pile-up fea- 
ture which results in a prolonged power-law up to the electron cut- 
off moment um p e ~ p max , where a steeper su per-exponential cutoff 
takes over (Zirakashvili & Aharonian 2007). Applying the theory 
of plane-parallel shock acceleration that is justified because of the 
large curvature radius of the shock, the equilibrium electron distri- 
bution function at each position in the shock is given by 



fc(X, Pe) = CePe [ 1 + j(.X, p c )f e eXp 



PmaxW 
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Here, x is a spatial coordinate along the shock normal, measured 
from the shock position and the electron cutoff momentum is 



PmaxW = (F + \X\GT 1 . 
_ U 



where 



and 
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, , (B8) 

The shape of the pile-up region of electrons in equation 1B5| | is 
given by S e = 9/5 and the power-law function 



4poaK 
ffinj + 2 a in j - 1 
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The upstream quantities have the index 1 and the downstream 
quantities the index 2. The flow velocity in the inertial frame of the 
shock is denoted by ;/ (u\ = v s , u 2 = ui/r c where r c = p 2 fp\ is 
the density compression factor at the shock), and p represents the 
injection momentum normalized with ra e c. Note that the electron 
cutoff momentum p max is independent of po, as expected. The ratio 
between the cooling and diffusive acceleration time scales given by 
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Here, the inverse energy loss time scale of an electron dominated 
by synchrotron and inverse-Compton losses is given by 



E 
E 



3m, c 



(s B + £ ph ) . 



(Bll) 



The CRe acceleration depends on the type of diffusion which we as- 
sume to be parallel to the magnetic field and to be in the Bohm-limit 
(as motivated by young supernova remnants observations, see also 
Section l2.4. It . The energy dependent diffusion constant is given by 



p c nteCr 



1- 
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which we assume to be the same across the shock, i.e. k — Ki - k 2 . 
We also use that the magnetic field fluctuations SB are of the same 
amplitude as the magnetic field rj = (B/SB) 2 ~ 1. 

Bohm diffusion becomes more effective at higher energies, 
which causes the cooling induced cutoff in the electron spectrum 
to move to lower energies as the electrons are transported advec- 
tively with the flow downstream. Integration over the post-shock 
volume causes the cutoffs to add up to the power-law that is steeper 
by unity compared to the injection power-law. Thus the steady-state 



spectrum that balance losses through IC/synchrotron cooling with 
gains through electron shock acceleration, is derived through the 
integral over the IC radiating volume Vol = SD (where S is the 
surface area, and D the thickness). It is given by 
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where erf is the standard Gaussian error function 
dAbramowitz & Stegunl 1 19651) . We do not expect the main 
characteristics of the pile-up region to change much when it is 
integr ated over the post-shock regime ( Zirakas hvili & Aharonian! 
2007); to obtain a consistent semi-analytic formula of the spec- 
trum, we determine the specific values of J(p c ), G 2 , and 5 e ,cooi 
through fits to numerically integrated spectra (see below). The 
distribution function in equation JB 1 3I > has a break in the spectrum 
at the electron momentum p e = Pe.break = 1/DG 2 and a cutoff at 
Pc = Pe.cut = 1 IF. The electron momentum cutoff is determined 
from the strength of the mag netic field a nd the properties of 
electron diffusion in the shock IWebb G.M.IIl984l) . and are given 
by 
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The equivalent electron energy cutoff in the relativistic regime is 
given by 

1O.5 
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In Fig. IB 1 1 we show a comparison between our electron distribu- 
ti on given by equatio n jB13b compared to the electron spectrum 
in lEnfilin et al.l l[l998) where an energy independent diffusion has 
been assumed. To plot the spectrum, we used typical values for the 
temperature, Mach number and density in the virial region. These 
values are found in Table [2] In addition we assumed a magnetic 
field amplification in the shock of a factor ten, resulting in a mag- 
netic field of about 6pG. The spectral comparison show that the 
spectral shapes are similar up to the cutoff region where the energy 
dependent Bohm diffusion in our model induces a steeper cutoff. 
Furthermore, the offset in the normalization of the two spectra be- 
tween the break and cutoff regime is caused by the integrated pile- 
up regime (equation[B9} that effectively increases the energy of the 
break from 1 /D G 2 to 1 /D G 2 ■ 

In Fig. IB2l we show both the cooled electron spectra that we 
derive by numerically integrating f c {x, p c ) over the shock as well 
as the fitted spectra. For / c (p e ) we use typical values for both 
the strength of the shocks (a c = 3.1) and the break in the elec- 
tron spectrum which we fix at a constant momentum (Pcbreak = 
300MeV/m e c 2 ). In our fit, we allow for three independent fit vari- 
ables, Ai, A 2 , and A3. We find that, 



A, = 1.45, A 3 = 1.95, where, 
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We note that the factor 4.5 in <5 e ,cooi follows from 
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Figure Bl. Comparison of the spectrum of primary electrons. We com- 
pare the electron spectrum in ou r model (equation |B13j using results from 
lzirakashvili_& Aharonian 1 2007) of the high-energy cutoff region to a spec- 
trum from lEnBlin et alj J 1998) where an energy independent diffusion has 
been assumed. There is a break in the spectrum at p c = 1 /DG2 and a cutoff 
at p c = 1 / F. The comparison shows similar spectral shape up to the cutoff 
region where the energy dependent Bohm diffusion in our model induces 
a much steeper cutoff. The offset of the two spectra between the break and 
cutoff regimes is caused by the integrated pile-up regime indicated by J(p c ) 
(equation IB9t that effectively shifts the break by a factor of two to higher 
energies. 



Zirakashvil i & Aharonianl J2007t) and represents the spectral 
flattening of the power-law in the low-energy regime up to the 
regime where the exponential cutoff dominates. 



B3 Inverse Compton radiation 

Inverse Compton scattering of CMB photons off ultra-relativistic 
electrons with Lorentz factors of y c > 10 4 redistributes these pho- 
tons into the hard X-ray/y-ray regime according to equation {6]l. 
The integrated IC source density A JC for an isotropic power-law 
distribution of CR electrons, as described by equation JB3b . is ob- 
ta ined by integrating th e IC so urce function s y {E y ) in equation (43) 
of IPfrommer & Enfilirj J2004» (in the case of Thomson scattering) 
over an energy interval between observed photon energies E\ and 
E 2 yielding 
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where a v = (a c - l)/2 denotes the photon spectral index, r c = 
e 2 /(m e c 2 ) the classical electron radius, and £(a) me Riemann f- 



£ m „ = 3 TeV 
£ m „= 10 TeV 
£„„ = 30 TeV 
R = 100 TeV 



£_ = 300 TeV 
£„„ = 1000 TeV 
No cutoff 



■S 10' r 




Figure B2. Cooled electron distribution of primary electrons for different 
cutoff energies. We show the electron distribution for a fixed spectral index 
a c = 3.1 and varying electron cutoff energies E nmK in raising order; 3 TeV 
(red), 10 TeV (yellow), 30 TeV (green), 100 TeV (light blue), 300 TeV 
(blue), 1000 TeV (purple), and no cutoff (black). The data (crosses) is gen- 
erated by numerically integrating equation )B13t for different cutoff ener- 
gies. The fit (solid lines) are derived using three free parameters, where the 
fitted variables are shown in equation IB 161 . 



function ( Abram owitz & St egun 1965). The CRe normalization C e 
is given by equation dB2t and dB4t for the secondary and primary 
electrons, respectively. 

In the following, we provide a simple analytic formula that 
captures the primary inverse-Compton emission from galaxy clus- 
ters for strong shocks and intermediate strength shocks (i.e. or in j ~ 
2 - 2.7) with an accuracy of five percent or better. In addition, we 
want the analytic formula to be valid in the full energy range of 
pIC emission, i.e. not limited by the Klein-Nishina (KN) suppres- 
sion where the center of mass energy of photons becomes compa- 
rable to the electron mass and the less efficient energy transfer from 
electron to photon causes a break in the photon spectrum. 

Using the exact spectra in the asymptotic low- and high- 
energy regime, together with the result of numerical calculations at 
intermediate energies, we parametrize the integrated source func- 
tion for pIC emission by, 
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where Bohm diffusion has been assumed. The normalization con- 
stants /to(£e,max> C e ) an d /ic(*e) ar e derived in equations dB201 l 
and l |B19l l, respectively. The shape of the IC spectrum without 
an exponential cutoff scales as E^ y in the Thomson regime, 
and ste epens to E^ c log(Ei C ) in th e KN suppressed high-energy 
regime dBlumenthal & Gould [l97(jl) . In Fig. IB3I we fit the numer- 
ically calculated intermediate energy regime using the following 
parametrization, 
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Figure B3. The primary inverse Compton y-ray number flux weighted by 
photon energy for different electron spectral indices and electron cutoff 
energies. Main panel: The data is numerically calculated using a cutoff en- 
ergy i? max = 3 TeV (crosses) and without a cutoff (X) for different electron 
spectral indices; a c = 2.6 (blue), o- c = 3.1 (green), a c = 3.6 (yellow), and 
a e =4.1 (red). The fits to the data using £ max = 3 TeV (no cutoff) are 
represented by the solid (dash-dotted) lines. Bottom panel: The difference 
between the data and the fits. The flux from the fits using the dominating 
electron spectral indices a c =3.1 and a e = 3.6, agree within a few percent 
with the data in the dominating flux regime. The flux from the fits using 
electron spectral indices of a e = 2.6 and a c =4.1 have slightly lower pre- 
cision, and agree within 10 percent with the data in the dominating flux 
regime. 
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which respects this asymptotic behavior of the IC spectrum. We 
use three free independent fit variables, £kn,i > jSkn, and Akn and 
find that 
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The resulting spectra accurately describe the pIC emission without 
an exponential cutoff for a e ~ 2.5 - 4.5. In addition, we capture the 
shape of the integrated source function for pIC in equation JB2U . 
where we include the super-exponential cutoff. We find that the 
shape of the transition region is well approximated by 



6ic(Eic,a e ) = 0.529tt e - 0.134 log,, 
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which depends on both the photon energy and the electron spectral 
index. In the process of finding a good fit, we allowed for multiple 
free parameters in equation dB21t . In the end, however, we only 
allowed die to vary and fixed all the other free parameters at typical 
values to keep the formula as simple as possible. 

In Fig. lB4l we test our analytic formula for pIC emission, given 
by equations dB 19t . JB20| > and JB21IB241 . to the numerically cal- 
culated pIC for different electron cutoff energies. We find that the 



Figure B4. The primary inverse Compton y-ray number flux weighted 
by photon energy for different electron cutoff energies. Main panel: The 
pIC emission are shown for a fixed electron spectral index (a e = 3.1) 
and varying electron cutoff energies (£ ma x) in raising order; 3 TeV (red), 
10 TeV (orange), 30 TeV (yellow), 100 TeV (green), 300 TeV (light blue), 
1000 TeV (blue), and no cutoff (black). The numerically calculated data 
is represented by crosses and the fits are shown with solid lines. Bottom 
panel: The difference between the data and the fits. The flux from the fits 
agree within a few percent with the data in the dominating flux regime. 



fits agree within a few percent with the numerically calculated IC 
emission in the dominating flux regime where the flux is larger than 
10 percent of the maximum pIC flux (that have been normalized 
with E'y'). 

This paper has been typeset from a TgX/ I5TjX file prepared by the 
author. 
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